Closed manishsoni0291 closed 1 year ago
barronh
commented
1 year ago This sounds like some sort of version issue. Run the command below, and send me the result.
python -c " import xarray as xr xr.show_versions() "
Also, please include the whole traceback. Right now, you are only showing the last 7 lines and the line that triggered it. There is a lot of code in between that I need to see to understand where this originates.
manishsoni0291
commented
1 year ago commit: None python: 3.9.6 (default, Jul 28 2021, 12:38:21) [GCC 4.8.5 20150623 (Red Hat 4.8.5-39)] python-bits: 64 OS: Linux OS-release: 3.10.0-1160.2.2.el7.x86_64 machine: x86_64 processor: x86_64 byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: ('en_US', 'UTF-8') libhdf5: 1.12.0 libnetcdf: 4.7.4
xarray: 2022.11.0 pandas: 1.3.1 numpy: 1.21.1 scipy: 1.7.0 netCDF4: 1.5.7 pydap: None h5netcdf: 1.0.2 h5py: 3.3.0 Nio: None zarr: None cftime: 1.5.0 nc_time_axis: None PseudoNetCDF: None rasterio: None cfgrib: None iris: None bottleneck: None dask: 2021.07.1 distributed: None matplotlib: 3.4.2 cartopy: None seaborn: None numbagg: None fsspec: 2021.07.0 cupy: None pint: None sparse: None flox: None numpy_groupies: 0.9.13 setuptools: 56.0.0 pip: 21.3.1 conda: None pytest: 6.2.4 IPython: 7.25.0 sphinx: 4.1.2
manishsoni0291
commented
1 year ago @barronh
barronh
commented
1 year ago I just pushed an update version 0.2.3 that updates the code as follows:
tm5_b = xr.DataArray(
tm5_constant_b.mean(1), dims=('layer',),
coords=[('layer', satl3f['layer'])]
)to
tm5_b = xr.DataArray(
tm5_constant_b.mean(1), dims=('layer',),
coords=dict(layer=satl3f['layer'])
)The dictionary form is not ambiguous, while the sequence of sequences is. I have verified that the old form raises an error with newer xarray versions (e.g., 0.20.0), and that the new code does not. I've also confirmed that it is backward compatible.
Can you confirm that the new version solves this problem?
manishsoni0291
commented
1 year ago @barronh The problem is solved. Now it is facing another error in terms of TSTEP and Date and time found. It is not able to find the TSTEP in netcdf file
header of TropOMINO2_2019-07-25_12US1_CMAQ.nc .............................................................................................................. netcdf TropOMINO2_2019-07-25_12US1_CMAQ { dimensions: LAY = 35 ; ROW = 299 ; COL = 459 ; variables: float NO2(LAY, ROW, COL) ; NO2:_FillValue = NaNf ; NO2:long_name = "NO2 " ; NO2:units = "ppmV " ; NO2:var_desc = "Variable NO2 " ; NO2:coordinates = "TSTEP" ;
................................................... header of ncdump -h CCTM_CONC_v52_noDDM_gcc_TRECH2_BOSTON_12km_base_20170725.nc
netcdf CCTM_CONC_v52_noDDM_gcc_TRECH2_BOSTON_12km_base_20170725 { dimensions: TSTEP = UNLIMITED ; // (25 currently) DATE-TIME = 2 ; LAY = 35 ; VAR = 179 ; ROW = 299 ; COL = 459 ; variables: int TFLAG(TSTEP, VAR, DATE-TIME) ; TFLAG:units = "<YYYYDDD,HHMMSS>" ; TFLAG:long_name = "TFLAG " ; TFLAG:var_desc = "Timestep-valid flags: (1) YYYYDDD or (2) HHMMSS " ; float NO2(TSTEP, LAY, ROW, COL) ; NO2:long_name = "NO2 " ; NO2:units = "ppmV " ; NO2:var_desc = "Variable NO2
Both are different.
barronh
commented
1 year ago I see that they are different, but that is different than there being an error. Is there an error?
manishsoni0291
commented
1 year ago Yep @barronh this is an error coz it doesn't find the TSTEP in the TropOMINO2_2019-07-25_12US1_CMAQ.nc i am supplying full header now. Can u check again the output by using your files.
barronh
commented
1 year ago this is an error coz it doesn't find the TSTEP
Just because it doesn't have a TSTEP does not mean there is an error. This file is not intended to be IOAPI compliant. So, it does not have a TSTEP dimension, nor does it have a TFLAG variable.
The file looks fine.
manishsoni0291
commented
1 year ago @barronh I tried opening the file in panoply and verdi it is giving error.
barronh
commented
1 year ago Manish,
Panoply support would be ideal, and I may look into improving compatibility with Panoply. Right now, Panoply is inferring that this is an IOAPI file and failing to read the time. The file is not intended to fully support IOAPI for many reasons. For example, the vertical structure is inherited from the satellite. The satellites can have vertical structures that are not supported by IOAPI.
So, I think your question is how can I make a figure? You can make a figure pretty easily without Panoply using xarray.
The commands below should make a nice little plot. I've chosen to visualize the CMAQ VCD and the TropOMI VCD processed with the CMAQ AMF. An alternative is to use the VCD from TropOMI and visualize CMAQ with the TropOMI averaging kernel. Either comparison is better than comparing CMAQ's VCD directly to TropOMI's.
import xarray as xr
import cmaqsatproc as csp
import matplotlib.pyplot as plt
# Open a GRIDDESC and a CMAQ processed TropOMI file
cg = csp.open_griddesc('12US1')
qf = xr.open_dataset('TropOMINO2_2019-07-25_12US1_CMAQ.nc')
# Make a figure to hold two panels.
fig, axx = plt.subplots(1, 2, figsize=(12, 4))
opts = dict(
norm=plt.Normalize(vmin=0),
cmap='viridis'
)
# Add CMAQ raw VCD
qf['VCDNO2_CMAQ'].plot(ax=axx[0], **opts)
# Add TropOMI VCD recalculated with CMAQ AMF
qf['VCDNO2_TOMI_CMAQ'].plot(ax=axx[1], **opts)
# add maps with states
for ax in axx.ravel():
cg.csp.cno.drawstates(ax=ax)
# Save the figure to disk
fig.savefig('vcd_no2.png')@barronh Hi Barron First of all thank for your assistance.
When i use the script to plot the figure. There is no difference coming. values are also very less.
When i use the same dimension plotting of figure completed successfully. When the CMAQ concentration file have lesser dimension (spatial reference) It gives error. it is shown below
import numpy
import pandas
import geopandas
import netCDF4
import h5netcdf
import pyproj
import cmaqsatproc as csp
import xarray as xr
GDNAM = '12NE3'
date='2018-07-01'
readername = 'TropOMINO2'
satreader = csp.reader_dict[readername]
l3 = xr.open_dataset(f'{readername}_{date}_{GDNAM}.nc')
qf = csp.open_ioapi(f'CCTM_CONC_4.454_Bench_2018_12NE3_20180701.nc')[['NO2']]
mf = csp.open_ioapi(f'MET_CRO_3D_4.454_Bench_2018_12NE3_20180701.nc')
qf['DENS'] = mf['DENS']
qf['ZF'] = mf['ZF']
qf['PRES'] = mf['PRES']
# Create satellite according to CMAQ, and CMAQ according to satellite
overf = satreader.cmaq_process(qf, l3)
ValueError Traceback (most recent call last)in 22 qf['PRES'] = mf['PRES'] 23 # Create satellite according to CMAQ, and CMAQ according to satellite ---> 24 overf = satreader.cmaq_process(qf, l3) ~/cmaqsatproc/cmaqsatproc/readers/tropomi/__init__.py in cmaq_process(cls, qf, satl3f, key) 678 overf['NO2_PER_M2'].attrs['units'] = 'mole/m**2'.ljust(16) 679 --> 680 ak = overf['NO2_AK_CMAQ'] = cls.cmaq_ak(overf, satl3f) 681 overf['NO2_SW_CMAQ'] = cls.cmaq_sw(overf, satl3f) 682 # uses AK for tropopause ~/cmaqsatproc/cmaqsatproc/readers/tropomi/__init__.py in cmaq_ak(cls, overf, satl3f, amfkey) 635 Averaging kernel on the CMAQ grid 636 """ --> 637 q_sw_trop = cls.cmaq_sw(overf, satl3f, amfkey=amfkey) 638 q_ak = q_sw_trop / satl3f['air_mass_factor_troposphere'] 639 q_ak.attrs.update(q_sw_trop.attrs) ~/cmaqsatproc/cmaqsatproc/readers/tropomi/__init__.py in cmaq_sw(cls, overf, satl3f, amfkey, tropopausekey) 578 Tropospheric scattering Weights on the CMAQ grid 579 """ --> 580 return TropOMI.cmaq_sw( 581 overf, satl3f, amfkey=amfkey, tropopausekey=tropopausekey 582 ) ~/cmaqsatproc/cmaqsatproc/readers/tropomi/__init__.py in cmaq_sw(cls, overf, satl3f, amfkey, tropopausekey) 339 sw_trop = ak * satl3f[amfkey] 340 --> 341 q_sw_trop = coord_interp( 342 overf['PRES'], pres, sw_trop, 343 indim='layer', outdim='LAY', ascending=False ~/cmaqsatproc/cmaqsatproc/utils.py in coord_interp(coordout, coordin, varin, verbose, interp, outdim, indim, ascending, **kwds) 276 XP = coordin.sortby(indim, ascending=ascending) 277 FP = varin.sortby(indim, ascending=ascending) --> 278 interped = xr.apply_ufunc( 279 interp1d, 280 FP, ~/.local/lib/python3.9/site-packages/xarray/core/computation.py in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args) 1202 # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc 1203 elif any(isinstance(a, DataArray) for a in args): -> 1204 return apply_dataarray_vfunc( 1205 variables_vfunc, 1206 *args, ~/.local/lib/python3.9/site-packages/xarray/core/computation.py in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args) 297 298 if len(args) > 1: --> 299 args = deep_align( 300 args, join=join, copy=False, exclude=exclude_dims, raise_on_invalid=False 301 ) ~/.local/lib/python3.9/site-packages/xarray/core/alignment.py in deep_align(objects, join, copy, indexes, exclude, raise_on_invalid, fill_value) 833 out.append(variables) 834 --> 835 aligned = align( 836 *targets, 837 join=join, ~/.local/lib/python3.9/site-packages/xarray/core/alignment.py in align(join, copy, indexes, exclude, fill_value, *objects) 770 fill_value=fill_value, 771 ) --> 772 aligner.align() 773 return aligner.results 774 ~/.local/lib/python3.9/site-packages/xarray/core/alignment.py in align(self) 557 self.find_matching_unindexed_dims() 558 self.assert_no_index_conflict() --> 559 self.align_indexes() 560 self.assert_unindexed_dim_sizes_equal() 561 ~/.local/lib/python3.9/site-packages/xarray/core/alignment.py in align_indexes(self) 402 if need_reindex: 403 if self.join == "exact": --> 404 raise ValueError( 405 "cannot align objects with join='exact' where " 406 "index/labels/sizes are not equal along " ValueError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'ROW' ('ROW',) ```
barronh
commented
1 year ago When the CMAQ concentration file have lesser dimension (spatial reference) It gives error. it is shown below
I do not understand this. What do you mean when the CMAQ concentration have lesser dimension?
The TropOMI file you make for 12NE3 must also be 12NE3 and the CONC file must be 3 dimensional.
If this does not clarify your question, can you post your files some where so that I can test it?
manishsoni0291
commented
1 year ago @barronh The problem was i haven't included the GRIDDESC in ~/cmaqsatproc/cmaqsatproc/cmaq.py
I also think this can also be included into documentation file for other users to tell the users to specify the GRIDDESC in cmaq.py. The comparison is not good in terms there is large underestimation. This is CMAQ benchmark case output taken from CMAQ site only using CMAQ 5.4. I have included the nc file as file as png file.
TropOMINO2_2018-07-01_2018_12NE3_CMAQ.nc.gz
barronh
commented
1 year ago I also think this can also be included into documentation file for other users to tell the users to specify the GRIDDESC in cmaq.py.
You don't need to include the GRIDDESC in the code. csp.open_griddesc has an optional keyword gdpath where you can point to your own custom file. I have, just now, improved the documentation on open_griddesc and other functions in cmaqsatproc.cmaq to make this obvious.
gf = csp.open_griddesc('12NE3', gdpath='/path/to/your/GRIDDESC')I think I know why your totals are low. It looks like your CONC file is just 1 layer tall. As a result, your vertical column density is only integrating one layer. You would need to re-run CMAQ with the full vertical structure.
p.s., I am going to make some edits to make this thread shorter and more readable for others in the future.
barronh
commented
1 year ago I am going to close this issue, since it was originally related to the xarray version compatibility. If you have more questions, please open a new issue.
I am getting this error when trying to invoke command
overf = satreader.cmaq_process(qf, l3)
elif isinstance(obj, tuple): 125 if isinstance(obj[1], DataArray): --> 126 raise TypeError( 127 "Using a DataArray object to construct a variable is" 128 " ambiguous, please extract the data using the .data property."
TypeError: Using a DataArray object to construct a variable is ambiguous, please extract the data using the .data property.