using CDO to extract dataset only for a specific region

Question

I want to use 'cdo' to extract data from a precipitation NetCDF dataset using another NetCDF over South America. I have tried multiple procedures, but I always get some error (such as grid size not same, Unsupported generic co-ordinates, etc).

The codes I have tried:

cdo mul chirps_2000-2015_annual_SA.nc Extract feature.nc output.nc
# Got grid error

cdo -f nc4 setctomiss,0 -gtc,0 -remapcon,r1440x720 Chirps_2000-2015_annual_SA.nc CHIRPS_era5_pev_2000-2015_annual_SA_masked.nc
# Got unsupported generic error

If the grid is missing then you might find data slicing in NCO more helpful, the links to the files are now dead so I removed the last sentence referring to the dead link. — ClimateUnboxed, Feb 17 '20 at 08:06

score 1 · Accepted Answer · answered Jun 18 '19 at 14:52

I am sure you could make/find more elegant solution, but I just combined Python and shell executable cdo to fulfill the task (calling subprocess can be considered as a bad habit sometimes/somewhere).

#!/usr/bin/env ipython
import numpy as np
from netCDF4 import Dataset
import subprocess
# -------------------------------------------------
def nc_varget(filein,varname):
    ncin=Dataset(filein);
    vardata=ncin.variables[varname][:];
    ncin.close()
    return vardata
# -------------------------------------------------
gridfile='extract_feature.nc'
inputfile='precipitation_2000-2015_annual_SA.nc'
outputfile='selected_region.nc'
# -------------------------------------------------
# Detect the start/end based on gridfile:
poutlon=nc_varget(gridfile,'lon')
poutlat=nc_varget(gridfile,'lat')

pinlon=nc_varget(inputfile,'lon')
pinlat=nc_varget(inputfile,'lat')

kkx=np.where((pinlon>=np.min(poutlon)) & (pinlon<=np.max(poutlon)))
kky=np.where((pinlat>=np.min(poutlat)) & (pinlat<=np.max(poutlat)))
# -------------------------------------------------
# -------------------------------------------------
commandstr='cdo selindexbox,'+str(np.min(kkx))+','+str(np.max(kkx))+','+str(np.min(kky))+','+str(np.max(kky))+' '+inputfile+' '+outputfile
subprocess.call(commandstr,shell=True)

The problem in your data is that the file "precipitation_2000-2015_annual_SA.nc" does not specify the grid at the moment - variables lon, lat are generic and hence the grid is generic. Otherwise you could use other operators instead of selindexbox. File extract_feature.nc are more closer to the standard as the variables lon, lat have also the name and unit attributes.

using CDO to extract dataset only for a specific region

1 Answers1