benbovy / PyGChem

A Python interface to the GEOS-Chem Model, a global 3-D chemical transport model (CTM) for atmospheric composition
GNU General Public License v3.0
10 stars 4 forks source link

Error when loading nested data from bpch files (iris backend) #1

Closed benbovy closed 8 years ago

benbovy commented 8 years ago

Loading nested data from bpch (using the iris backend) gives the following error:

ValueError                                Traceback (most recent call last)
<ipython-input-12-2cbf0abd887c> in <module>()
----> 1 datasets.load( a[0] )
/usr/local/anaconda-2.3.0/lib/python2.7/site-packages/iris/__init__.pyc in load(uris, constraints, callback)
    282 
    283     """
--> 284     return _load_collection(uris, constraints, callback).merged().cubes()
    285 
    286 
/usr/local/anaconda-2.3.0/lib/python2.7/site-packages/iris/__init__.pyc in _load_collection(uris,    constraints, callback)
    252     try:
    253         cubes = _generate_cubes(uris, callback, constraints)
--> 254         result = iris.cube._CubeFilterCollection.from_cubes(cubes, constraints)
    255     except EOFError as e:
    256         raise iris.exceptions.TranslationError(
/usr/local/anaconda-2.3.0/lib/python2.7/site-packages/iris/cube.pyc in from_cubes(cubes, constraints)
    134         pairs = [_CubeFilter(constraint) for constraint in constraints]
    135         collection = _CubeFilterCollection(pairs)
--> 136         for cube in cubes:
    137             collection.add_cube(cube)
    138         return collection
/usr/local/anaconda-2.3.0/lib/python2.7/site-packages/iris/__init__.pyc in _generate_cubes(uris, callback, constraints)
    239         if scheme == 'file':
    240             part_names = [x[1] for x in groups]
--> 241             for cube in iris.io.load_files(part_names, callback, constraints):
    242                 yield cube
    243         elif scheme in ['http', 'https']:
/usr/local/anaconda-2.3.0/lib/python2.7/site-packages/iris/io/__init__.pyc in load_files(filenames, callback, constraints)
    194                 yield cube
    195         else:
--> 196             for cube in handling_format_spec.handler(fnames, callback):
    197                 yield cube
    198 
/usr/local/anaconda-2.3.0/lib/python2.7/site-packages/PyGChem-0.3.0dev-py2.7.egg/pygchem/datafield_backends/iris_backend.pyc in cubes_from_bpch(filenames, callback, **kwargs)
     81             if leading_datablock['origin'] != (1, 1, 1):
     82                 imin = np.array(leading_datablock['origin']) - 1
---> 83                 imax = imin + np.array(leading_datablock['shape'])
     84                 region_box = zip(imin, imax)
     85             else:
ValueError: operands could not be broadcast together with shapes (3,) (2,) 
tsherwen commented 8 years ago

This issue appears to only be linked to 0.5_0.666 resolution. I now have GEOS-Chem nested running at 0.25_0.3125 and PyGchem/Iris has no issues reading outputted ctm.bpch files. Issue with 0.5*0.666 ctm.bpch files remains apparent. 11137139_10153075270726434_1256669371465821_n

benbovy commented 8 years ago

Actually, the issue seems to appear when loading a bpch file in which the 1st datablock is a 2-dimensional field with origin != (1, 1, 1).

When loading a bpch file into iris cubes, the 1st datablock in the file was taken as reference in order to avoid creating new (duplicate) iris coordinates for each datablock.

It should now be fixed in fa0d2e59ea4d666a416751ba16fbcf110c8294f5. I implemented a basic cache system to create and assign coordinates to new created cubes.

Can you tell me if it is working with your data (use the pygchem's master branch) before I close the issue?

tsherwen commented 8 years ago

Thanks for taking a look.

On updating to the latest PyGChem master branch. The previous error message no longer appears, however a new error is present ( see output below).

The script being used to test this: https://github.com/tsherwen/AC_tools/blob/master/bpch2netCDF.py

Traceback (most recent call last): File "species_surface_mean_val_3D.py", line 21, in debug=debug ) File "/Users/Tomas/Google Drive/PhD/PhD_progs/AC_tools/funcs4GEOSC.py", line 1294, in get_GC_output convert_to_netCDF( wd ) File "/Users/Tomas/Google Drive/PhD/PhD_progs/AC_tools/bpch2netCDF.py", line 25, in convert_to_netCDF data = datasets.load(bpch_files) File "//anaconda/lib/python2.7/site-packages/iris/init.py", line 284, in load return _load_collection(uris, constraints, callback).merged().cubes() File "//anaconda/lib/python2.7/site-packages/iris/init.py", line 254, in _load_collection result = iris.cube._CubeFilterCollection.from_cubes(cubes, constraints) File "//anaconda/lib/python2.7/site-packages/iris/cube.py", line 136, in from_cubes for cube in cubes: File "//anaconda/lib/python2.7/site-packages/iris/init.py", line 241, in _generate_cubes for cube in iris.io.load_files(part_names, callback, constraints): File "//anaconda/lib/python2.7/site-packages/iris/io/init.py", line 196, in load_files for cube in handling_format_spec.handler(fnames, callback): File "//anaconda/lib/python2.7/site-packages/pygchem/dataset_backends/backend_iris.py", line 115, in cubes_from_bpch dim_coords = _get_datablock_dim_coords(datablock, coord_cache) File "//anaconda/lib/python2.7/site-packages/pygchem/dataset_backends/backend_iris.py", line 63, in _get_datablock_dim_coords og, sp = leading_datablock['origin'], leading_datablock['shape'] NameError: global name 'leading_datablock' is not defined

benbovy commented 8 years ago

Oups, fixed typo in ba8976bab9bfd7b9f61bf3c80f7cefb8bf48307c. Sorry!

tsherwen commented 8 years ago

Great. Everything seems to be working for 0.5*0.666 now. Thanks Ben.