geoschem / gcpy

Python toolkit for GEOS-Chem. Contains basic plotting scripts, plus the suite of GEOS-Chem benchmarking utilities.
https://gcpy.readthedocs.io
Other
50 stars 24 forks source link

[BUG/ISSUE] Issue in gcpy.compare_single_level: Top title is not printed #1

Closed yantosca closed 5 years ago

yantosca commented 5 years ago

I was using gcpy.compare_single_level to compare output between GCC and GCHP runs. The data variable I am using is Met_MODISLAI (which is only xy dimension). The data variables are shown here:

Ref (GCC)
<xarray.Dataset>
Dimensions:       (lat: 91, lon: 144, time: 1)
Coordinates:
  * time          (time) datetime64[ns] 2016-01-15T01:00:00
  * lon           (lon) float64 -180.0 -177.5 -175.0 ... 172.5 175.0 177.5
  * lat           (lat) float64 -89.5 -88.0 -86.0 -84.0 ... 84.0 86.0 88.0 89.5
Data variables:
    Met_MODISLAI  (time, lat, lon) float32 ...

Dev
<xarray.Dataset>
Dimensions:       (lat: 288, lon: 48, time: 1)
Coordinates:
  * lon           (lon) float64 1.0 2.0 3.0 4.0 5.0 ... 44.0 45.0 46.0 47.0 48.0
  * lat           (lat) float64 1.0 2.0 3.0 4.0 5.0 ... 285.0 286.0 287.0 288.0
  * time          (time) datetime64[ns] 2016-01-15T01:00:00
Data variables:
    Met_MODISLAI  (time, lat, lon) float32 ...

And my calling sequence is:

#!/usr/bin/env python

# Imports
import os
import xarray as xr
import gcpy

# Ref dataset
refstr  = '12.3.0 GCC'
refdir  = '~/LT/12.3.0/geosfp_2x25_TransportTracers/OutputDir'
reffile = 'GEOSChem.StateMet_inst.20160115_0100z.nc4'
refpath = os.path.join(refdir, reffile)
refdata = xr.open_dataset(refpath)

# Dev dataset
devstr  = '12.3.0 GCHP'
devdir  = '~/SP/if17/gchp_TransportTracers/OutputDir'
devfile = 'GCHP.StateMet_inst.20160115_0100z.nc4'
devpath = os.path.join(devdir, devfile)
devdata = xr.open_dataset(devpath)

varlist = ['Met_MODISLAI']

# Variables to plot

# PDF filename
pdfname = '/n/scratchlfs/jacob_lab/ryantosca/GC/rundirs/plots/pdf/MODIS_LAI_surface_nc_diag.pdf'

# Compare data
gcpy.compare_single_level(refdata, refstr, devdata, devstr,
                          varlist=varlist, pdfname=pdfname)

It produces the PDF, but for some reason the top title is not printed. The echoback to stdout is:

Reuse existing file: ./conservative_2x2.5_1x1.25.nc
/net/seasasfs02/srv/export/seasasfs02/share_root/ryantosca/python/geo/miniconda/envs/geo/ ib/python3.6/site-packages/xesmf/backend.py:36: UserWarning: Input array is not F_CONTIGUUS. Will affect performance.
  warnings.warn("Input array is not F_CONTIGUOUS. "
Reuse existing file: ./conservative_c48_1x1.25_0.nc
Reuse existing file: ./conservative_c48_1x1.25_1.nc
Reuse existing file: ./conservative_c48_1x1.25_2.nc
Reuse existing file: ./conservative_c48_1x1.25_3.nc
Reuse existing file: ./conservative_c48_1x1.25_4.nc
Reuse existing file: ./conservative_c48_1x1.25_5.nc

Creating /n/scratchlfs/jacob_lab/ryantosca/GC/rundirs/plots/pdf/MODIS_LAI_surface_nc_diagpdf for 1 variables
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-ryantosca'
/net/seasasfs02/srv/export/seasasfs02/share_root/ryantosca/python/geo/miniconda/envs/geo/ib/python3.6/site-packages/xarray/core/dataarray.py:507: FutureWarning: elementwise compaison failed; returning scalar instead, but in the future will perform elementwise comparion
  return key in self.data
0 Incorrect dimensions for Met_MODISLAI!

See attached file below.

I believe that the top title is computed in this block:

        if 'lev' in refdata[varname].dims and 'lev' in devdata[varname].dims:
            if ilev == 0: levstr = 'Surface'
            elif ilev == 22: levstr = '500 hPa'
            else: levstr = 'Level ' +  str(ilev-1)
            figs.suptitle('{}, {}'.format(varname,levstr), fontsize=fontsize, y=offset)
        elif 'lat' in refdata[varname].dims and 'lat' in devdata[varname] and 'lon' in refdata[varname].dims and 'lon' in devdata[varname]: 
            figs.suptitle('{}'.format(varname), fontsize=fontsize, y=offset)
        else:
            print('Incorrect dimensions for {}!'.format(varname))   

So I am not sure what is going on here.

no_toptitle

yantosca commented 5 years ago

This issue was caused by a typo. The line:

        elif 'lat' in refdata[varname].dims and 'lat' in devdata[varname] and 'lon' in   refdata[varname].dims and 'lon' in devdata[varname]:

should have been:

         elif 'lat' in refdata[varname].dims and 'lat' in devdata[varname]**.dims** and 'lon' in   refdata[varname].dims and 'lon' in devdata[varname]**.dims**:

This is now fixed in commit #0d60595.