2

I am trying to extract data from each grid cell in a user-specified (non rectangular) domain defined by lat / lon boundaries. My input file is on a curvilinear grid. I have tried various methods python ,cdo, ncks, but I still cannot figure it out. I just want a timeseries of information for each grid cell in my polygon domain subset of the input ncfile. My input ncfile information is given here:

$ ncdump -h 1979_sfc_out.nc 

netcdf \1979_sfc_out {
dimensions:
    x = 83 ;
    y = 94 ;
    time = UNLIMITED ; // (8736 currently)
    nv4 = 4 ;
variables:
    float time(time) ;
        time:axis = "T" ;
        time:long_name = "time" ;
        time:standard_name = "time" ;
        time:units = "hours since 1979-1-2 00:00:00" ;
        time:calendar = "standard" ;
    float x(x) ;
        x:axis = "x" ;
        x:long_name = "X-coordinate in Cartesian system" ;
        x:standard_name = "projection_x_coordinate" ;
        x:units = "meters" ;
    float y(y) ;
        y:axis = "y" ;
        y:long_name = "Y-coordinate in Cartesian system" ;
        y:standard_name = "projection_y_coordinate" ;
        y:units = "meters" ;
    float lon(y, x) ;
        lon:units = "degrees_east" ;
        lon:valid_range = -180., 180. ;
        lon:standard_name = "longitude" ;
        lon:bounds = "lon_bnds" ;
    float lat(y, x) ;
        lat:units = "degrees_north" ;
        lat:valid_range = -90., 90. ;
        lat:standard_name = "latitude" ;
        lat:bounds = "lat_bnds" ;
    float lon_bnds(y, x, nv4) ;
        lon_bnds:units = "degreesE" ;
    float lat_bnds(y, x, nv4) ;
        lat_bnds:units = "degreesN" ;
    char mapping ;
        mapping:false_easting = 0. ;
        mapping:false_northing = 0. ;
        mapping:grid_mapping_name = "polar_stereographic" ;
        mapping:latitude_of_projection_origin = 90. ;
        mapping:standard_parallel = 64. ;
        mapping:straight_vertical_longitude_from_pole = -152. ;
        mapping:semi_major_axis = 6370000. ;
        mapping:semi_minor_axis = 6370000. ;
    float SEAICE(time, y, x) ;
        SEAICE:_FillValue = -9999.f ;
        SEAICE:units = "fraction" ;
        SEAICE:long_name = "Ice concentration (ice=1;no ice=0)" ;
        SEAICE:grid_mapping = "mapping" ;
        SEAICE:coordinates = "lon lat" ;

Some of the things I have tried are

lat1=71.2
lat2=72.9
lon1=-176.5
lon2=-160

cdo sellonlatbox,lon1,lon2,lat1,lat2 $ifile $box1_ofile
cdo sellonlatbox (Abort): Float parameter >lon1< contains invalid character at position 1!

I think the problem is my input file has the x,y dimension in meters (which does not have a negative sign and likely is the 'invalid character at position 1'), and I am asking cdo to extract the lat/lon dimension in degrees. I have the variable 'mapping' in my ncfile which might be useful to convert meters to lat/lon, but I cannot figure out how to do it.

And ncks doesn't work for me here, likely due to the same meter <-> lon/lat problem I am seeing with cdo.

ncks -v SEAICE,U10,V10 -d latitude,71.2,72.9 -d longitude,-176.5,-160. $ifile -O $box1_ofile

ncks: ERROR dimension latitude is not in input file

Although I want a polygon extracted, these examples I have tried are just rectangle subsets, which I thought to get the polygon I can do multiple rectangle subsets to achieve my final polygon shape, but if there is a better way to do this any advice is appreciated.

Thanks

ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86
pwprnt
  • 521
  • 3
  • 9
  • 27

1 Answers1

2

NCO's auxiliary coordinates feature are designed to hyperslab curvilinear coordinates on an unstructured grid with 1-D lat and lon. Try this:

ncks -X lon_min,lon_max,lat_min,lat_max in.nc out.nc
ncks -X -176.5,-160.,71.2,72.9 in.nc out.nc

You may also daisy-chain multiple -X options to get weird polygons. Unfortunately that probably won't work for you because, I just noticed, you have a curvilinear grid with 2D lat and lon. For that try the ncap2 where feature. The examples in the manual show how to mask along rectangular boundaries, and you can daisy chain the conditions in the where() statement to obtain a polygon. This will not change the dimensionality of the output file, but it allows you to set everything outside the polygon to a _FillValue.

Charlie Zender
  • 5,929
  • 14
  • 19