0

Given two nc files here:

data/CMIP6_UKESM1-0-LL_Lmon_piControl_r1i1p1f2_gpp_1960-3059.nc
data/CMIP6_UKESM1-0-LL_Amon_piControl_r1i1p1f2_tas_1960-3059.nc

Read the first file:

from netCDF4 import Dataset
import numpy as np

ds1 = Dataset('data/CMIP6_UKESM1-0-LL_Lmon_piControl_r1i1p1f2_gpp_1960-3059.nc')
print(ds1.variables.keys()) # get all variable names

Out:

odict_keys(['gpp', 'time', 'time_bnds', 'lat', 'lat_bnds', 'lon', 'lon_bnds', 'clim_season', 'season_year'])

Read the second file:

ds2 = Dataset('data/CMIP6_UKESM1-0-LL_Amon_piControl_r1i1p1f2_tas_1960-3059.nc')
print(ds2.variables.keys()) 

Out:

odict_keys(['tas', 'time', 'time_bnds', 'lat', 'lat_bnds', 'lon', 'lon_bnds', 'clim_season', 'height', 'season_year'])

Check gpp variable:

gpp = ds1.variables['gpp'] # gpp variable
print(gpp)

Out:

<class 'netCDF4._netCDF4.Variable'>
float32 gpp(time, lat, lon)
    _FillValue: 1e+20
    standard_name: gross_primary_productivity_of_biomass_expressed_as_carbon
    long_name: Carbon Mass Flux out of Atmosphere Due to Gross Primary Production on Land [kgC m-2 s-1]
    units: kg m-2 s-1
    cell_methods: area: mean where land time: mean
    coordinates: clim_season season_year
unlimited dimensions: 
current shape = (3300, 144, 192)
filling on

Check tas variable:

tas = ds2.variables['tas'] # tas variable
print(tas)

Out:

<class 'netCDF4._netCDF4.Variable'>
float32 tas(time, lat, lon)
    _FillValue: 1e+20
    standard_name: air_temperature
    long_name: Near-Surface Air Temperature
    units: K
    cell_methods: area: time: mean
    coordinates: clim_season height season_year
unlimited dimensions: 
current shape = (3300, 144, 192)
filling on

Now I would like to calculate correlation between gpp and tas, then plot their correlation values on the map.

How could I do that? Thanks.

ah bon
  • 9,293
  • 12
  • 65
  • 148

1 Answers1

1

You should be able to do this easily using my package nctoolkit (https://nctoolkit.readthedocs.io/en/latest/).

My understanding is that you want to plot the temporal correlation coefficient per grid-cell. In that case:

import nctoolkit as nc
files = ["data/CMIP6_UKESM1-0-LL_Lmon_piControl_r1i1p1f2_gpp_1960-3059.nc",
"data/CMIP6_UKESM1-0-LL_Amon_piControl_r1i1p1f2_tas_1960-3059.nc"]
ds = nc.open_data(files)
ds.merge()
ds.cor_time(var1 = "gpp", var2 =  "tas")
ds.plot()

If you wanted the spatial correlation coefficient between variables per time step, you can use the following line in place of the second last one:

ds.cor_space(var1 = "gpp", var2 =  "tas")
Robert Wilson
  • 3,192
  • 11
  • 19
  • Sorry, it raises an error: `AttributeError: 'LGEOS360' object has no attribute 'GEOSBufferWithParams'`. – ah bon Apr 15 '21 at 09:12
  • this sounds like a cartopy issue with conda https://stackoverflow.com/questions/65027227/attributeerror-lgeos360-object-has-no-attribute-geosbufferwithparams – Robert Wilson Apr 15 '21 at 09:48
  • 1
    `ds.to_xarray().cor.plot()` should work in place of the final line for a quick plot, to get around any cartopy issues – Robert Wilson Apr 15 '21 at 10:08
  • The error comes at first line of code: `nc.open_data`. – ah bon Apr 15 '21 at 10:09
  • 1
    Do you mean the `import nctoolkit` line? There is nothing in `open_data` that should throw an error like that. I was working with UKESM cmip6 data just last week, so the package should be able to handle it problem free – Robert Wilson Apr 15 '21 at 11:08
  • It raise a new error: `ValueError: The files in the dataset to do not have the same grid. Consider using regrid!` Did you test with my data from the link? – ah bon Apr 15 '21 at 15:26
  • Yes. I've tested it problem free. Which version of Python are using? I only recently realized conda-forge has stopped building nctoolkit for 3.8 for macos. The package itself is tested daily on macos with 3.8, so problem free. If you are on 3.8 you should maybe do a pip install – Robert Wilson Apr 15 '21 at 15:56
  • Could you please upload the output map in your answer? – ah bon Apr 16 '21 at 05:45
  • One more question, if we have more than two files, how could we calculate their correlation two by two? – ah bon Apr 16 '21 at 05:53
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/231200/discussion-between-robert-wilson-and-ah-bon). – Robert Wilson Apr 16 '21 at 08:18