barronh / pseudonetcdf

PseudoNetCDF like NetCDF except for many scientific format backends
GNU Lesser General Public License v3.0
76 stars 35 forks source link

CAMx tutorial #39

Closed MPeasley42 closed 6 years ago

MPeasley42 commented 6 years ago

Hi. I've been tasked to teach my modeling group to use CAMx. I am attempting to run the TCEQ simulation , but am getting this error message ' File "marsha_test.py", line 30 print "Gridded emission file variable keys", gefile.variables.keys() ^ SyntaxError: invalid syntax '

Since I've only just learned RedHat Linux, I'm at a complete loss as to what needs to be changed in the python script. Any help will be much appreciated.

MPeasley42 commented 6 years ago

Hi @barronh Scratch the previous post. I've solved that one. This is my new problem : Traceback (most recent call last): File "marsha_test.py", line 18, in windfile = wind(windpath, rows = 65, cols = 83) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/camxfiles/wind/Memmap.py", line 70, in init self.__time_hdr_fmts={12: "fii", 8: "fi"}[rf.record_size] KeyError: 622940243

barronh commented 6 years ago

Hi @MPeasley42

Thanks for your interest in PseudoNetCDF. You request for help provides little information that would allow me to respond. For example, there is no context. What version of CAMx inputs are you trying to read, what is the domain etc. Which version of PseudoNetCDF did you install, etc... I'll do my best.

I suspect you're running a relatively recent TCEQ simulation. That tells me two things. (1) you are likely using the wrong reader and (2) the rows and columns you are using are from an old TCEQ domain -- none of the new domains are 65x83.[1]

As of CAMx v6, the wind reader is no longer needed. Wind inputs, as of v6, are written in the more self-describing uamiv format. Thus, you should import the uamiv reader, which does not require that you know the dimensions of the data. The code might now look like

# or from PseudoNetCDF.camxfiles.Memmaps import uaimv
from PseudoNetCDF import pnc
windfile = pnc('--format=uamiv', windpath)

Followed by whatever else you want to do.

[1] https://www.tceq.texas.gov/airquality/airmod/data/domain

MPeasley42 commented 6 years ago

Hi @barronh I will try my best to give as much info as I can. I've been kind of thrown into this and I'm learning as I go. We are using CAMx v6. I am using the TCEQ simulation in the CAMx tutorial. As for which reader I am to use or domains I am to have, I'm very foggy on that. I've just been given a server, a crash course on RedHat, and now, I have to sink or swim so to speak. Thank you for your input. I will keep posted as I make the changes.

MPeasley42 commented 6 years ago

Hi @barronh I'm at a complete lost and have no idea what I'm doing wrong this is my python script :

import libraries

from PseudoNetCDF.camxfiles.Memmaps import from numpy import

or from PseudoNetCDF.camxfiles.Memmaps import uaimv

from PseudoNetCDF import pnc

set file paths

windpath = 'camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.v45' zppath = 'camx_zp.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.v45' kvpath = 'camx_kv.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.TKE.v45' lupath = 'camx_landuse.hgbpa_04km.utcsr' temppath = 'camx_temp.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.v45' crpath = 'camx_cr.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.v45' humpath = 'camx_hum.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.v45' ptpath = 'camx_cb05_ei_el.20050519.hgb8h2.bc05may.reg7a_pscfv2' gepath = 'camx_cb05_ei_lo.20050519.hgb8h2.bc05may.reg8.hgbpa_04km' aconcpath = 'camx451_cb05.20050519.hgb8h2.bc05may.reg8_pscfv2.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.avrg.hgbpa_04km'

Open files

windfile = pnc('--format=uamiv', windpath) zpfile = height_pressure(zppath, rows = 191, cols = 218) kvfile = vertical_diffusivity(kvpath, rows = 191, cols = 218) lufile = landuse(lupath, rows = 191, cols = 218) tempfile = temperature(temppath, rows = 191, cols = 218) crfile = cloud_rain(crpath, rows = 191, cols = 218) humfile = humidity(humpath, rows = 191, cols = 218) ptfile = point_source(ptpath) gefile = uamiv(gepath) aconcfile = uamiv(aconcpath)

Perform demonstration calculations

print - "Gridded emission file variable keys", gefile.variables.keys() print - "Gridded NO emission rate unit", gefile.variables['NO'].units print - "Gridded emission NOx average over space in %s" % gefile.variables['NO'].units, eval('NO[:] + NO2[:]', globals(), gefile.variables).sum(3).sum(2).sum(1) print - "Maximum pressure in %s" % zpfile.variables['PRES'].units, zpfile.variables['PRES'][:].max() print - "Average layer 1, noon temperature in %s" % tempfile.variables['AIRTEMP'].units, tempfile.variables['AIRTEMP'][12, 0].mean() print - "Maximum cloud water content in %s" % crfile.variables['CLOUD'].units, crfile.variables['CLOUD'][:].mean() print - "Median wind speed in %s" % windfile.variables['U'].units, median(eval('sqrt(U2 + V2)', globals(), windfile.variables)) print - "Ozone Changes larger than 20ppb/hr", eval('diff(O3, axis = 0) > 0.02', globals(), aconcfile.variables).sum()

R = 831.4472 # m3 hPa K-1 mol-1 Na = 6.02214179e23 # avagadro's number molecules/mol dZ = eval('concatenate([HGHT[:,[0]], diff(HGHT[:], axis = 1)], axis = 1)', globals(), zpfile.variables) # m vol = dZ gefile.XCELL gefile.YCELL # m3 temp = tempfile.variables['AIRTEMP'][:] # K pres = zpfile.variables['PRES'][:] # hPa print - "Min air number density (molecules/m*3)", (pres vol / temp / R * Na).min()

--- END EXAMPLE

And this is the resulting response in my terminal: Traceback (most recent call last): File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/pncparse.py", line 665, in getfiles f = allreaders[file_format](ipath, **format_options) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/camxfiles/uamiv/Memmap.py", line 98, in init self.readheader() File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/camxfiles/uamiv/Memmap.py", line 199, in readheader GDTYPE = self.GDTYP={0: 1, 1: 5, 2: 2, 3: 6}[cproj] KeyError: 1147496033

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "marsha_test.py", line 23, in windfile = pnc('--format=uamiv', windpath) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/pncparse.py", line 572, in pnc out = PNC(*args, ifiles = ifiles, actions = actions, kwds).ifiles File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/pncparse.py", line 557, in PNC ifiles, outargs = pncparse(args = args, ifiles = ifiles, parser = parser, kwds) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/pncparse.py", line 419, in pncparse ifiles.extend(pncprep(subarg)[0]) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/pncparse.py", line 611, in pncprep fs = getfiles(ipaths, args) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/pncparse.py", line 670, in getfiles raise oute# from e OSError: Unable to open path with uamiv(path, **{}) path="camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_fddats_uhsst_utcsrlulc.v45" error="1147496033" As I said, I've been thrown into this with no real training and I haven't a clue as how to proceed.

barronh commented 6 years ago

@MPeasley42 -

I did a little google. It looks like you are trying to follow a blog post of mine from 2011. Since that data is from the 2005 episode with the old domains, the number of rows and columns was actually right in your first message. Since the 2005 episode used CAMx v4.5, the reader was also right.

I downloaded the file an example file from TCEQ and the code below works just fine. A few tips, you are using python3, so print statements are now print functions (e.g., print 1 is now print(1)).

import numpy as np
from PseudoNetCDF.camxfiles.Memmaps import wind
windfile = wind('camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45', rows = 65, cols = 83)
print("Median wind speed in %s" % windfile.variables['U'].units, float(np.median(eval('np.sqrt(U**2 + V**2)', globals(), windfile.variables))))

If that doesn't work with your wind file, I suggest you redownload it.

Let's take a huge step backward.

1) Whose research group are you in? University? PI? 2) What exactly have you been tasked with? 2a) Are you tasked with running a contemporary CAMx run or TCEQ's 2005 simulation? 2b) Are you tasked with teaching the group to run CAMx or analyze inputs?

-Barron

MPeasley42 commented 6 years ago

Hi @barronh I'm part of the Air Dispersion Modeling group with the Division for Air Quality in the state of Kentucky. Previously, we depended on EPA region 4 in Atlanta, Georgia for our regional haze and ozone modeling. Over the years, we have found that our actual impacts were lower than Atlanta's numbers due do our particular weather patterns and terrain. Also, our need for more realistic SO2 and NOx impacts than what was given from AERMOD called for the used of photochemical modeling. Since I am the only one in my group with some programming back ground, I am tasked to learn how to run CAMx along with all of the necessary pre and post processors. Afterwards, I am to train the rest of my group. And I thought being a new mom was stressful 😄.

-Marsha

MPeasley42 commented 6 years ago

Hello @barronh I've updated the script and my response was : Traceback (most recent call last): File "marsha_test.py", line 26, in windfile = wind('camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45', rows = 65, cols = 83) File "/home/mpeasley/anaconda3/lib/python3.6/site-packages/PseudoNetCDF/camxfiles/wind/Memmap.py", line 70, in init self.__time_hdr_fmts={12: "fii", 8: "fi"}[rf.record_size] KeyError: 622940243 I'm not sure what this error is. I've found the file and line 70. I'm not sure of what to change to fix it.

Marsha

barronh commented 6 years ago

@MPeasley42 -

There is not an error in the reader. I think the problem is the file. In fact, I get a very similar error if I try to read the wind file without unzipping it.

Below are three code snippets to help you.

The first line calls md5sum and the second confirms the output. I suspect your md5sum doesn't match.

$ md5sum camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45
$ stat camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45

If your file is good, you should get exactly these values. c94182d214bfb821a7b65a5c94e46369 camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45 and a size from stat of 30224000

The script below uses Python3 to download an unzip the file.

from urllib.request import urlretrieve
import gzip
import shutil

urlwind = 'ftp://amdaftp.tceq.texas.gov/HGB8H2/camx/basecase/bc05may.reg10.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell/input/met/camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45.gz'
gzpath = 'camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45.gz'
windpath = 'camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45'
gzf = urlretrieve(urlwind, gzpath)
with gzip.open(gzpath, 'rb') as f_in:
    with open(windpath, 'wb') as f_out:
        shutil.copyfileobj(f_in, f_out)

alternatively, you can use a shell script

$ wget ftp://amdaftp.tceq.texas.gov/HGB8H2/camx/basecase/bc05may.reg10.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell/input/met/camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45.gz
$ gunzip camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45.gz

Once your md5sum matches exactly, you can run the script below.

from PseudoNetCDF.camxfiles.Memmaps import wind
import numpy as np
windpath = 'camx_wind.20050519.hgbpa_04km.2005ep0_eta_dbemis_upsfcfddats_newuhsst_newutcsrlulc_grell.v45'
windfile = wind(windpath, rows = 65, cols = 83)
print("Median wind speed in %s" % windfile.variables['U'].units, float(np.median(eval('np.sqrt(U**2 + V**2)', globals(), windfile.variables))))

I have run these scripts and have no problems. I'll need you to take it from here for a bit.

Given your level of experience, I'd recommend that you take a CAMx class and use VERDI or Panoply for visualization to get started. VERDI and Panoply are interactive point and click tools. Neither will let you read an old wind file, but no one should be using those any more anyway.

-B

MPeasley42 commented 6 years ago

Hi @barronh After much discussion with my management, they have decided to go with VERDI. Thank you so much for your help. _Marsha