Scaling a Subset of Input Data in NetCDF Format

2021-02-10·
Ka Ming FUNG
Ka Ming FUNG
· 11 min read

Reference: For general python handling of nc files using Xarray: http://xarray.pydata.org/en/stable/index.html

I am working on a set of input data of Ocean DMS concentration. I suspect the reported values over the Southern Ocean (vaguely defined as ocean south to 30ºS) is a bit high. I want to reduce those values by 50% to see if my model will match observations better.

Surprisingly, I didn’t find a quick way to do it. Here I present a work-a-round. (Please let me know if you have a smarter way.)

import wget                        # for downloading data via URL links
import xarray as xr                   # the major tool to work with NetCDF data!
from matplotlib import pyplot as plt  # for plotting

Open the monthly Ocean DMS nc file

The marine DMS data used in this example can be downloaded here. Or, we can download it from Python.

wget.download('http://rain.ucis.dal.ca/ctm/HEMCO/DMS/v2015-07/DMS_lana.geos.1x1.nc')

'DMS_lana.geos.1x1.nc'

We will use the function xr.open_dataset from the xarray package.

ds = xr.open_dataset('DMS_lana.geos.1x1.nc')

ds # display the dataset

<xarray.Dataset>
Dimensions:    (lat: 181, lon: 360, time: 12)
Coordinates:


    

  • lon        (lon) float32 -180.0 -179.0 -178.0 -177.0 … 177.0 178.0 179.0
    • 

  • lat        (lat) float32 -89.75 -89.0 -88.0 -87.0 … 87.0 88.0 89.0 89.75
    • 

  • time       (time) datetime64[ns] 1985-01-01 1985-02-01 … 1985-12-01
Data variables:
DMS_OCEAN  (time, lat, lon) float32 …
Attributes:
CDI:          Climate Data Interface version 1.6.4 (http://code.zmaw.de/p…
Conventions:  COARDS
history:      Mon Jul 20 16:29:35 2015: ncatted -O -a units,DMS_OCEAN,o,c…
Title:        COARDS/netCDF file created by BPCH2COARDS (GAMAP v2-17+)
Format:       NetCDF-3
Model:        GEOS5_47L
Delta_Lon:    1.0
Delta_Lat:    1.0
NLayers:      47
Start_Date:   19850101
Start_Time:   0
End_Date:     19850201
End_Time:     0
Delta_Time:   744
CDO:          Climate Data Operators version 1.6.4 (http://code.zmaw.de/p…

We can also visualize annual mean of the data for a quick check:

ds['DMS_OCEAN'].mean(dim=['time']).plot()

<matplotlib.collections.QuadMesh at 0x7fd61bbd01f0>

png

Scaling down the Southern Ocean DMS

As mentioned, we want to scale down the ocean DMS below 30ºS.

# Step 1: make two copies of the data set, one for grid cells with lat > 30ºS, one for <30ºS

dms_above30S = ds['DMS_OCEAN'].sel(lat=slice(-30,90))

dms_below30S = ds['DMS_OCEAN'].sel(lat=slice(-90,-30-1e-9)) # we add a small value here to avoid the grid cells with lat = 30ºS to appear in both data sets.

# Step 2: reduce the DMS in the dataset of interest by 50% 

dms_below30S = dms_below30S * 0.5

# Step 3: merge the untouched and scaled data
ds_partly_scaled = xr.merge([dms_below30S, dms_above30S], compat="no_conflicts")

To ensure our exercise is correct, let’s make some plots to compare the DMS concentrations in ds_partly_scaled and the original dataset, ds.

(ds_partly_scaled['DMS_OCEAN'] / ds['DMS_OCEAN']).plot(col='time', col_wrap=3)

<xarray.plot.facetgrid.FacetGrid at 0x7fd61bd3ba00>

png

Yellow = 1 (untouched); Purple = 0.5(scaled). Mission accomplished?

ds_partly_scaled # missing attributes from the original dataset.

<xarray.Dataset>
Dimensions:    (lat: 181, lon: 360, time: 12)
Coordinates:


    

  • 

    lat        (lat) float64 -89.75 -89.0 -88.0 -87.0 … 87.0 88.0 89.0 89.75
    

    • 

  • 

    lon        (lon) float32 -180.0 -179.0 -178.0 -177.0 … 177.0 178.0 179.0
    

    • 

  • 

    time       (time) datetime64[ns] 1985-01-01 1985-02-01 … 1985-12-01
Data variables:
DMS_OCEAN  (time, lat, lon) float32 0.0 0.0 … 1.4103641e-07 1.4103641e-07

Preserving metadata from the original nc file

Not yet! Now the new merged data set is missing the Attributes of the original.

To fix this, we need to make a xarray copy of the original to preserve the original information, but with the new scaled values.

ds_new = ds.copy(data=ds_partly_scaled)

ds_new # now with the Attributes

<xarray.Dataset>
Dimensions:    (lat: 181, lon: 360, time: 12)
Coordinates:


    

  • 

    lon        (lon) float32 -180.0 -179.0 -178.0 -177.0 … 177.0 178.0 179.0
    

    • 

  • 

    lat        (lat) float64 -89.75 -89.0 -88.0 -87.0 … 87.0 88.0 89.0 89.75
    

    • 

  • 

    time       (time) datetime64[ns] 1985-01-01 1985-02-01 … 1985-12-01
Data variables:
DMS_OCEAN  (time, lat, lon) float32 0.0 0.0 … 1.4103641e-07 1.4103641e-07
Attributes:
CDI:          Climate Data Interface version 1.6.4 (http://code.zmaw.de/p…
Conventions:  COARDS
history:      Mon Jul 20 16:29:35 2015: ncatted -O -a units,DMS_OCEAN,o,c…
Title:        COARDS/netCDF file created by BPCH2COARDS (GAMAP v2-17+)
Format:       NetCDF-3
Model:        GEOS5_47L
Delta_Lon:    1.0
Delta_Lat:    1.0
NLayers:      47
Start_Date:   19850101
Start_Time:   0
End_Date:     19850201
End_Time:     0
Delta_Time:   744
CDO:          Climate Data Operators version 1.6.4 (http://code.zmaw.de/p…

Finally, save the scaled DMS in a new file

And, it’s one-line easy.

ds_new.to_netcdf('DMS_lana.geos.1x1.SouthernHalved.nc') # always good to have a descriptive name

Again, I know that it’s likely not the best way. Please let me know if you have a smarter way. I am always open to learn!

Python
Ka Ming FUNG
Authors
Ka Ming FUNG
Lead Data Scientist