GCPy not properly recognizing 0.125x0.15625 lat-lon grid

[msulprizio]

Name and Institution (Required)

Name: Melissa Sulprizio
Institution: Harvard

Description of your issue or question

When trying to create plots comparing 12-km (0.125x0.15625 degree) lat-lon grid to other resolutions, I am getting a similar error as reported in #313.

ValueError: operands could not be broadcast together with shapes (402,447) (402,416)

I am able to plot the 0.125x0.15625 output against itself, but I noticed the resolution string is incorrect. For example:

I added debug statements in grid.py for the longitude:
[-129.984 -129.8277 -129.6714... -60.5868 -60.4305 -60.2742]

which returns the longitude resolution of
0.15629999999998745

I suspect there's an issue with numerical roundoff happening and we haven't encountered it because previous resolution strings haven't extended beyond 4 decimal places.

[yantosca]

@msulprizio: I will look into this. It may be an easy fix like you said, due to roundoff.

[yantosca]

@msulprizio: I believe the issue is in the get_input_res function of gcpy/grid.py. This is used to get the lon & lat resolution from the input xr.Dataset or xr.DataArray to the plotting codes.

It could be that the roundoff is happening in the netCDF files themselves. I tested the algorithm in get_input_res with this script:

#!/usr/bin/env python3

import numpy as np
from gcpy.grid import make_grid_LL

grid_12km = make_grid_LL("0.125x0.15625")

grid_12km["lat"].sort()
grid_12km["lon"].sort()
lat_res = np.abs(grid_12km["lat"][2] - grid_12km["lat"][1])
lon_res = np.abs(grid_12km["lon"][2] - grid_12km["lon"][1])

res = str(lat_res) + "x" + str(lon_res)

print(lon_res)
print(lat_res)
print(res)

and I got this output

0.15625
0.125
0.125x0.15625

So I think the algorithm is OK. Maybe the the longitude and latitude values in the netCDF file have some roundoff (if they are REAL*4) and this could be causing it. We could put in a workaround to force the correct resolution string.

[yantosca]

@msulprizio: Has this been resolved?

[msulprizio]

I have not confirmed that this is resolved now. I will report back when I get the 0.125x0.15625 transport tracer runs past a wetdep error.

[msulprizio]

This is still an issue. I am comparing 0.125x0.15625 output over North America to 0.25x0.3125 using the compare_diags.py. script and this compare_diags.yml as input. I am still getting the error:

Using configuration file compare_diags.yml                                                                                                                                                                                    

Comparing variable names in compare_varnames                                                                                                                                                                                  
32 common variables                                                                                                                                                     
0 variables in ref only                                                                                                                                                                 
0 variables in dev only                                                                                                                                                                 
All variables have same dimensions in ref and dev                                                                                                                                                                             
... Creating single-level plots                                                                                                                                                                                         
joblib.externals.loky.process_executor._RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/externals/loky/process_executor.py", line 463, in _process_worker
    r = call_item()
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/externals/loky/process_executor.py", line 291, in __call__
    return self.fn(*self.args, **self.kwargs)
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 589, in __call__
    return [func(*args, **kwargs)
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 589, in <listcomp>
    return [func(*args, **kwargs)
  File "/n/home05/msulprizio/python/gcpy/gcpy/plot/compare_single_level.py", line 794, in createfig
    absdiff = np.array(ds_dev_cmp) - np.array(ds_ref_cmp)
ValueError: operands could not be broadcast together with shapes (402,448) (401,448) 
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/n/holylfs05/LABS/jacob_lab/Users/msulprizio/Runs/12km_validation/test_1day/./compare_diags.py", line 367, in <module>
    main(sys.argv)
  File "/n/holylfs05/LABS/jacob_lab/Users/msulprizio/Runs/12km_validation/test_1day/./compare_diags.py", line 362, in main
    compare_data(config, read_data(config))
  File "/n/holylfs05/LABS/jacob_lab/Users/msulprizio/Runs/12km_validation/test_1day/./compare_diags.py", line 292, in compare_data
    compare_single_level(
  File "/n/home05/msulprizio/python/gcpy/gcpy/plot/compare_single_level.py", line 1150, in compare_single_level
    results = Parallel(n_jobs=n_job)(
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 1952, in __call__
    return output if self.return_generator else list(output)
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 1595, in _get_outputs
    yield from self._retrieve()
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 1699, in _retrieve
    self._raise_error_fast()
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 1734, in _raise_error_fast
    error_job.get_result(self.timeout)
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 736, in get_result
    return self._return_or_raise()
  File "/n/home05/msulprizio/python/micromamba/envs/gcpy_env/lib/python3.9/site-packages/joblib/parallel.py", line 754, in _return_or_raise
    raise self._result
ValueError: operands could not be broadcast together with shapes (402,448) (401,448)

[yantosca]

@msulprizio: I wonder if this is a halfpolar grid thing. I think the nested-grid output doesn't use halfpolar=true (or something like that).

[github-actions]

This issue has been automatically marked as stale because it has not had recent activity. If there are no updates within 7 days it will be closed. You can add the "never stale" tag to prevent the issue from closing this issue.