“”” .. CMAQ Spatial Subset
CMAQ Spatial Subset
This example shows how to spatially subset a file in two ways. The first subset windows a file and retains the source file resolution. The second subset resamples the coarse data using cell centroid to extract concentrations from the source file and put it in a finer resolution file.
Both examples use he extent of the 12US1 file defined in GRIDDESC to create a bounding box and apply it to the 108NHEMI2 domain data (CONC.nc).
CMAQ Window
This first code block takes a polar stereographic file and extracts just the cells within the 12US1 coninental domain. The meta-data is updated automatically, so this is pretty easy.
import PseudoNetCDF as pnc
# create a GRIDDESC
with open('GRIDDESC', 'w') as gf:
gf.write(
"' '\n'LamCon_40N_97W'\n 2 33.000 45.000 -97.000 -97.000 40.000\n" +
"' '\n'12US1'\n'LamCon_40N_97W' " +
"-2556000.0 -1728000.0 12000.0 12000.0 459 299 1\n' '"
)
# Create a Template
gf = pnc.pncopen('GRIDDESC', format='griddesc')
concf = pnc.pncopen('CONC.nc', format='ioapi').copy()
i, j = concf.ll2ij(
gf.variables['longitude'][:],
gf.variables['latitude'][:]
)
islice = slice(i.min(), i.max() + 1)
jslice = slice(j.min(), j.max() + 1)
wndwf = concf.sliceDimensions(COL=islice, ROW=jslice)
wndwf.save('wndw.nc', format='NETCDF3_CLASSIC')
CMAQ Resampling
This second code block is a continuation of the one above. Instead of windowing and preserving source resolution, we are resampling to produce a file of the target resolution. The extraction method is complex and so much of the code is spent making sure the meta-data is kept up to date.
# Get values for each i/j
tmpf = concf.sliceDimensions(LAY=0).sliceDimensions(COL=i, ROW=j, newdims=('ROW', 'COL'))
startdate = concf.getTimes()[0]
# Create an output file for resampled data
gf.updatetflag(startdate=startdate, overwrite=True)
gf.VGLVLS = tmpf.VGLVLS
nsteps = len(tmpf.dimensions['TSTEP'])
resamplef = gf.subsetVariables(['DUMMY']).sliceDimensions(
TSTEP=[0]*nsteps
)
# Set times to be sequential
resamplef.updatetflag(startdate=startdate, tstep=concf.TSTEP, overwrite=True)
# Copy in the data
for key, var in tmpf.variables.items():
if key == 'TFLAG': continue
resamplef.copyVariable(var, key=key)
del resamplef.variables['DUMMY']
delattr(resamplef, 'VAR-LIST')
resamplef.updatemeta()
# Make sure TFLAG matches
resamplef.updatetflag(tstep=concf.TSTEP, overwrite=True)
# Save file
resamplef.save('resample.nc')