PpCO2 for gn models

[abby-baskind]

Disclaimer: This is a very chaotic issue because getting the gn models to work is taking a backseat to making my presentation right now. I probably could have tried to debug better, but I mostly needed to post the issue before I forgot.

So I'm trying to plot meridional slices of PpCO2 for the gn models I regridded, like I did for the 'gr' models. This is the gr output I'm trying to replicate with gn.

I regridded the gn models like this...

targetgrid_ds = xe.util.grid_global(1.0, 1.0)

dd_regrid={}
for name,item in dd_new_new_gn.items():
    regridder = xe.Regridder(item, targetgrid_ds, 'bilinear', 
                         periodic=True, ignore_degenerate=True)
    ds_regridded = regridder(item)
    dd_regrid[name]=ds_regridded
list(dd_regrid.keys())

And I calculated PpCO2 with these functions.

def calc_PpCO2(ds):
    p = gsw.p_from_z(-1*ds['lev'], ds['y'], geo_strf_dyn_height=0, sea_surface_geopotential=0)
    insitutemp = gsw.t_from_CT(ds['so'], ds['thetao'], p)  # 0 is just a filler rn 
    conversion =  1e6/1035
    results = pyco2.sys(par1=ds.talk*conversion,par2=ds.dissic*conversion,par1_type=1,par2_type=2,
                        pressure_out=0, temperature_out = ds.thetao, pressure = p, 
                        temperature = insitutemp)
    return ds['talk'].copy(data=results['pCO2_out'])

def meridionalsection(ax,da,clims=None,title=None):
    im = ax.pcolormesh(da['y'],da['lev'],da)
    if clims is not None:
        im.set_clim(clims)
#     ax.set_xlim([-80,60])
    ax.set_xlim([20,160])
    ax.invert_yaxis()
    plt.colorbar(im,ax=ax)
    ax.set_title(title)

def calc_sigma2(ds):
    return gsw.sigma2(ds['so'],ds['thetao'])

def meridionalsection_with_sigma2(ax,da,sigma2,clims=None,title=None):
    meridionalsection(ax,da,clims)
    ax.contour(da['y'],da['lev'],sigma2,levels=[36,36.4,36.8],colors='w')
    ax.set_title(title)

I tried a few different ways. First, I tried plotting the models after they had been regridded like this

fig_pco2, axarr_pco2 = plt.subplots(nrows = 5, ncols = 1, figsize=[15,12])
fig_pco2.tight_layout(pad = 3.5)
plt.rc('font', size = 12)
plt.rc('axes', titlesize= 10)    
plt.rc('axes', labelsize= 10)
plt.rc('figure', titlesize=14)

ax_idx = 0 # index for your axes array
for name, ds_pco2 in dd_regrid.items():
    ax = axarr_pco2.flat[ax_idx]
    ds_pco2 = ds_pco2.isel(time=0).sel(x=slice(180,200)).mean('x',keep_attrs=True)
    pco2 = calc_PpCO2(ds_pco2)
    sigma2 = calc_sigma2(ds_pco2)
    meridionalsection_with_sigma2(ax,pco2,sigma2, clims=[500,2000],title=name)
    ax_idx += 1 

And got this ugly output.

I thought maybe the x slice I used was wrong, since longitude might not be 0 to 360 and might be -180 to 180 instead. So I changed slice(180, 200) to slice(-160, -180). Admittedly, I might have converted the longitudes wrong.

fig_pco2, axarr_pco2 = plt.subplots(nrows = 5, ncols = 1, figsize=[15,12])
fig_pco2.tight_layout(pad = 3.5)
plt.rc('font', size = 12)
plt.rc('axes', titlesize= 10)    
plt.rc('axes', labelsize= 10)
plt.rc('figure', titlesize=14)

ax_idx = 0 # index for your axes array
for name, ds_pco2 in dd_regrid.items():
    ax = axarr_pco2.flat[ax_idx]
    ds_pco2 = ds_pco2.isel(time=0).sel(x=slice(-160,-180)).mean('x',keep_attrs=True)
    pco2 = calc_PpCO2(ds_pco2)
    sigma2 = calc_sigma2(ds_pco2)
    meridionalsection_with_sigma2(ax,pco2,sigma2, clims=[500,2000],title=name)
    ax_idx += 1 

But the output was even worse

so then I tried plotting PpCO2 without regridding

fig_pco2, axarr_pco2 = plt.subplots(nrows = 5, ncols = 1, figsize=[15,12])
fig_pco2.tight_layout(pad = 3.5)
plt.rc('font', size = 12)
plt.rc('axes', titlesize= 10)    
plt.rc('axes', labelsize= 10)
plt.rc('figure', titlesize=14)

ax_idx = 0 # index for your axes array
for name, ds_pco2 in dd_new_new_gn.items():
    ax = axarr_pco2.flat[ax_idx]
    ds_pco2 = ds_pco2.isel(time=0).sel(x=slice(180,200)).mean('x',keep_attrs=True)
    pco2 = calc_PpCO2(ds_pco2)
    sigma2 = calc_sigma2(ds_pco2)
    meridionalsection_with_sigma2(ax,pco2,sigma2, clims=[500,2000],title=name)
    ax_idx += 1 

The output was questionable

But then I tried regridding after I calculated PpCo2

targetgrid_ds = xe.util.grid_global(1.0, 1.0)
# regridder = xe.Regridder(dd_new_new, targetgrid_ds, 'bilinear', 
#                          periodic=True, ignore_degenerate=True)

regridder = xe.Regridder(item, targetgrid_ds, 'bilinear', 
                         periodic=True, ignore_degenerate=True)

fig_pco2, axarr_pco2 = plt.subplots(nrows = 5, ncols = 1, figsize=[15,12])
fig_pco2.tight_layout(pad = 3.5)
plt.rc('font', size = 12)
plt.rc('axes', titlesize= 10)    
plt.rc('axes', labelsize= 10)
plt.rc('figure', titlesize=14)

ax_idx = 0 # index for your axes array
for name, ds_pco2 in dd_new_new_gn.items():
    ax = axarr_pco2.flat[ax_idx]
    ds_pco2 = ds_pco2.isel(time=0).sel(x=slice(180,200)).mean('x',keep_attrs=True)
    pco2 = calc_PpCO2(ds_pco2)
    sigma2 = calc_sigma2(ds_pco2)
    pco2_regrid = regridder(pco2)
    sigma2_regrid = regridder(sigma2)
    meridionalsection_with_sigma2(ax,pco2_regrid,sigma2_regrid, clims=[500,2000],title=name)
    ax_idx += 1 

And got this error

-------------
AssertionError                            Traceback (most recent call last)
<ipython-input-6-f73155172a99> in <module>
     22     pco2 = calc_PpCO2(ds_pco2)
     23     sigma2 = calc_sigma2(ds_pco2)
---> 24     pco2_regrid = regridder(pco2)
     25     sigma2_regrid = regridder(sigma2)
     26     meridionalsection_with_sigma2(ax,pco2_regrid,sigma2_regrid, clims=[500,2000],title=name)

/srv/conda/envs/notebook/lib/python3.8/site-packages/xesmf/frontend.py in __call__(self, indata, keep_attrs)
    414             return self.regrid_dask(indata)
    415         elif isinstance(indata, xr.DataArray):
--> 416             return self.regrid_dataarray(indata, keep_attrs=keep_attrs)
    417         elif isinstance(indata, xr.Dataset):
    418             return self.regrid_dataset(indata, keep_attrs=keep_attrs)

/srv/conda/envs/notebook/lib/python3.8/site-packages/xesmf/frontend.py in regrid_dataarray(self, dr_in, keep_attrs)
    464         input_horiz_dims, temp_horiz_dims = self._parse_xrinput(dr_in)
    465 
--> 466         dr_out = xr.apply_ufunc(
    467             self._regrid_array,
    468             dr_in,

/srv/conda/envs/notebook/lib/python3.8/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)
   1126     # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc
   1127     elif any(isinstance(a, DataArray) for a in args):
-> 1128         return apply_dataarray_vfunc(
   1129             variables_vfunc,
   1130             *args,

/srv/conda/envs/notebook/lib/python3.8/site-packages/xarray/core/computation.py in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args)
    269 
    270     data_vars = [getattr(a, "variable", a) for a in args]
--> 271     result_var = func(*data_vars)
    272 
    273     if signature.num_outputs > 1:

/srv/conda/envs/notebook/lib/python3.8/site-packages/xarray/core/computation.py in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)
    722             )
    723 
--> 724     result_data = func(*input_data)
    725 
    726     if signature.num_outputs == 1:

/srv/conda/envs/notebook/lib/python3.8/site-packages/xesmf/frontend.py in _regrid_array(indata, weights, shape_in, shape_out, sequence_in)
    426         if sequence_in:
    427             indata = np.expand_dims(indata, axis=-2)
--> 428         return apply_weights(weights, indata, shape_in, shape_out)
    429 
    430     @property

/srv/conda/envs/notebook/lib/python3.8/site-packages/xesmf/smm.py in apply_weights(weights, indata, shape_in, shape_out)
    102     extra_shape = indata.shape[0:-2]
    103 
--> 104     assert shape_horiz == shape_in, (
    105         'The horizontal shape of input data is {}, different from that of'
    106         'the regridder {}!'.format(shape_horiz, shape_in)

AssertionError: The horizontal shape of input data is (75, 332), different from that ofthe regridder (291, 360)!

And it makes sense that the shapes don't match since I'm averaging over x, but I don't know how to fix it. I'm also not fully convinced that changing the shape of the regridder or pco2 will fix the issue

[jbusecke]

The regridding afterwards will probably not work.

And got this ugly output.

Why ugly? Those look at least pretty consistent? You might end up in the wrong latitude given by how much land is in the North? To check this can you plot a map (of the first time step) of the regridded data? And a printout of the full regridded dataset with dimensions and coordinates?

I think you might need to assign some coordinates to the regridded data. if the values are just logical indicies (0-how long your x-dimension is) then the .sel command will not have the desired results. Fixing that should be relatively easy.

[gmacgilchrist]

I would agree that it's not necessarily ugly, just something going wrong in the sub selecting of data.

Question: does the regridding work as expected for T and S (and DIC and Alk)? I.e. is this an issue with the regridding, or is it somehow linked to the calculation of a secondary variable (PpCO2)? I would be quite surprised if it were the latter - the PpCO2 calculation should not be doing anything particularly special to the data structures.

If it's not related to the PpCO2 calculation, then this step is adding a layer of complexity that could make the regridding harder to debug.

[gmacgilchrist]

I thought maybe the x slice I used was wrong, since longitude might not be 0 to 360 and might be -180 to 180 instead. So I changed slice(180, 200) to slice(-160, -180). Admittedly, I might have converted the longitudes wrong.

You could check this directly by looking at the x coordinate.

Incidentally, you are selecting the longitudes here based on the x coordinate, as opposed to the 2D longitude array (i.e. what you had previously done with using the .where command). Is that because the regridded data is on a regular lat-lon grid?

[gmacgilchrist]

In the meridionalsection function you are selecting a latitude range (here specified as the xlims of the axis) between 20 and 160. Clearly these are not the actual latitudes of the data, but more likely the indices, as @jbusecke suggests.

.assign_coords applied to the dataset(s) or dataarrays will allow you to replace the indices with the actual latitude of the data. Then you should be able to select the appropriate latitudes for plotting.

This is likely also the case for the longitude dimension (the y dimension).