[GIS] netCDF-java has trouble reading this netCDF file(format) and fetching values from the latlon dimension

javanetcdfraster

I'm developing an application where I need to extract variable-values from netCDF files.

I'm processing several different formatted netCDF files and now one format seems to be creating some issues.

The problem is to correctly fetch data from lat/lon variables. Some lat/lon works, but the mapping does not seem to be complete.

ncdump of the netCDF file:

netcdf test2 {
dimensions:
        x = 182 ;
        y = 149 ;
        lev = 31 ;
        time = UNLIMITED ; // (72 currently)
variables:
        double lev(lev) ;
                lev:long_name = "generic" ;
                lev:units = "level" ;
                lev:axis = "Z" ;
        double time(time) ;
                time:standard_name = "time" ;
                time:units = "months since 1991-01-15 12:00:00" ;
                time:calendar = "standard" ;
        float variable1(time, lev, y, x) ;
                variable1:long_name = "Variable 1" ;
                variable1:units = "mg" ;
                variable1:_FillValue = 9.96921e+036f ;
                variable1:coordinates = "nav_lon nav_lat time lev" ;
        float nav_lon(y, x) ;
                nav_lon:long_name = "longitude" ;
                nav_lon:units = "degrees_east" ;
        float nav_lat(y, x) ;
                nav_lat:long_name = "latitude" ;
                nav_lat:units = "degrees_north" ;

// global attributes:
                :CDI = "Climate Data Interface version 1.5.1 (http://code.zmaw.de/projects/cdi)" ;
                :Conventions = "CF-1.0" ;
                :Title = "Test" ;
                :CDO = "Climate Data Operators version 1.5.1 (http://code.zmaw.de/projects/cdo)" ;
data:

 lev = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
    21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31 ;

 time = 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121,
    122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135,
    136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
    150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163,
    164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177,
    178, 179 ;
}

I can also show you subset of the nav_lon data (from opendap) if that could be of help to understand the format:

Dataset {
    Float32 nav_lon[y = 4][x = 4];
} test2.nc;
---------------------------------------------
nav_lon[4][4]
[0], 78.00004, 80.0, 81.99996, 83.99992
[1], 78.00004, 80.0, 81.99996, 83.99992
[2], 78.00004, 80.0, 81.99996, 83.99992
[3], 35.53008, 36.53008, 37.53008, 38.53008

Now to the actual code that is not producing the desired results. I'm using the latest version of NetCDF-Java, 4.3.19:

public class TestReadNetCDF {
    public static void main(String[] args) {
        String fileURI = "URI_TO_THE_FILE";
        readNetCDFFile(fileURI);
    }

    private static void readNetCDFFile(String filename) {
        NetcdfDataset ncfile = null;
        try {
            ncfile = NetcdfDataset.openDataset(filename);
            GridDataset grid = new GridDataset(ncfile);
            System.out.println("GRID:" + grid);
            process(grid);
        } catch (IOException ioe) {
            log("trying to open " + filename, ioe);
        } finally {
            if (null != ncfile)
                try {
                    ncfile.close();
                } catch (IOException ioe) {
                    log("trying to close " + filename, ioe);
                }
        }
        NetcdfDataset.shutdown();
    }
    private static void process(GridDataset gds) throws IOException {
        System.out.println("VARS:");
        for (VariableSimpleIF v : gds.getDataVariables())
            System.out.println(v.getFullName());
        GridDatatype grid = gds.findGridDatatype("variable1");
        GridCoordSystem gcs = grid.getCoordinateSystem();
        double lat = -63.63;
        double lon = -14.83;
        //lat = -61.24769;
        //lon = 78.00004;
        lat = -77.7742;
        lon = 78.00004;

        // find the x,y index for a specific lat/lon position
        // xy[0] = x, xy[1] = y
        int[] xy = gcs.findXYindexFromLatLonBounded(lat, lon, null);

        // read the data at that lat, lon and the first time and z level (if any)
        // note order is t, z, y, x
        System.out.println("x,y:" + xy[0] + "," + xy[1]);
        Array data = grid.readDataSlice(0, 0, xy[1], xy[0]);
        // we know its a scalar
        double val = data.getDouble(0);

        System.out.printf("Value at %f %f == %f%n", lat, lon, val);
    }
}

The problem seems to be that findXYindexFromLatLonBounded gets the wrong x/y indices (on most input values). If I drop "bounded", then it will fetch -1, -1, even though the latlon are within the "range" – though the mapping does not seem to be complete.

I could mention that ncWMS can read these files just fine and I can display the files as a WMS (and get point data out from there, but that is really inefficient and I'd like something better). Thredds can also read the files and output them through i.e. OpenDAP or Godiva2. If I use the OpenDAP URL it provides the same result as if I use the nc-file URI.

Any thoughts?

Best Answer

I cant tell if the lat / lon 2D arrays are correct, or if you are seeing a bug.

One way to test is to bring it up in the ToolsUI (example given).

Go to FeatureType / Grids, open the dataset and see what the grid looks like by hitting the red alien button (twice).

Publishing the dods URL lets others look at it remotely; you should also send what version of netcdf-java you are using.

Related Solutions

Python – Subsetting a Curvilinear netCDF File Using a lon/lat Bounding Box

If you want to subset some data from a NetCDF file using a lon/lat bounding box and that NetCDF file is not aligned with east/north, one strategy is to use a point-in-polygon routine and then find the min/max i,j indices of those points to define a subset to extract.

There are several different python packages to deal specifically with netCDF data from ROMS (a few examples are: https://github.com/kshedstrom/pyroms, and https://github.com/hetland/octant) but as far as I know they don't include subsetting by longitude/latitude. So here's an example:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import path 
import netCDF4

def bbox2ij(lon,lat,bbox=[-160., -155., 18., 23.]):
    """Return indices for i,j that will completely cover the specified bounding box.     
    i0,i1,j0,j1 = bbox2ij(lon,lat,bbox)
    lon,lat = 2D arrays that are the target of the subset
    bbox = list containing the bounding box: [lon_min, lon_max, lat_min, lat_max]

    Example
    -------  
    >>> i0,i1,j0,j1 = bbox2ij(lon_rho,[-71, -63., 39., 46])
    >>> h_subset = nc.variables['h'][j0:j1,i0:i1]       
    """
    bbox=np.array(bbox)
    mypath=np.array([bbox[[0,1,1,0]],bbox[[2,2,3,3]]]).T
    p = path.Path(mypath)
    points = np.vstack((lon.flatten(),lat.flatten())).T   
    n,m = np.shape(lon)
    inside = p.contains_points(points).reshape((n,m))
    ii,jj = np.meshgrid(xrange(m),xrange(n))
    return min(ii[inside]),max(ii[inside]),min(jj[inside]),max(jj[inside])

url = 'http://geoport.whoi.edu/thredds/dodsC/coawst_4/use/fmrc/coawst_4_use_best.ncd'
nc = netCDF4.Dataset(url)
#list variables
nc.variables.keys()
lon_rho = nc.variables['lon_rho'][:]
lat_rho = nc.variables['lat_rho'][:]
bbox = [-71., -63.0, 41., 44.]
i0,i1,j0,j1 = bbox2ij(lon_rho,lat_rho,bbox)
h = nc.variables['h'][j0:j1, i0:i1]
lon = lon_rho[j0:j1, i0:i1]
lat = lat_rho[j0:j1, i0:i1]
land_mask = 1 - nc.variables['mask_rho'][j0:j1, i0:i1]

# plot up subset
plt.pcolormesh(lon,lat,np.ma.masked_where(land_mask,-h),vmin=-350,vmax=0)
plt.colorbar()
plt.show()

which produces this plot here:

enter image description here

Update: I fixed the namespace issues reported by Luke below (I knew I would get burned by my ipython notebook --pylab eventually), and it occurs to me that I should probably post an example of subsampling velocity in ROMS, because it's trickier. Here I'm using the shrink and rot2d functions from Rob Hetland's Octant package to average the velocities to grid cell centers and rotate for plotting. Basemap isn't necessary here, but that's the example I have lying around, so here it is:

def rot2d(x, y, ang):
    '''rotate vectors by geometric angle

    This routine is part of Rob Hetland's OCTANT package:
    https://github.com/hetland/octant
    '''
    xr = x*np.cos(ang) - y*np.sin(ang)
    yr = x*np.sin(ang) + y*np.cos(ang)
    return xr, yr

def shrink(a,b):
    """Return array shrunk to fit a specified shape by triming or averaging.

    a = shrink(array, shape)

    array is an numpy ndarray, and shape is a tuple (e.g., from
    array.shape). a is the input array shrunk such that its maximum
    dimensions are given by shape. If shape has more dimensions than
    array, the last dimensions of shape are fit.

    as, bs = shrink(a, b)

    If the second argument is also an array, both a and b are shrunk to
    the dimensions of each other. The input arrays must have the same
    number of dimensions, and the resulting arrays will have the same
    shape.

    This routine is part of Rob Hetland's OCTANT package:
    https://github.com/hetland/octant

    Example
    -------

    >>> shrink(rand(10, 10), (5, 9, 18)).shape
    (9, 10)
    >>> map(shape, shrink(rand(10, 10, 10), rand(5, 9, 18)))        
    [(5, 9, 10), (5, 9, 10)]   

    """

    if isinstance(b, np.ndarray):
        if not len(a.shape) == len(b.shape):
            raise Exception, \
                  'input arrays must have the same number of dimensions'
        a = shrink(a,b.shape)
        b = shrink(b,a.shape)
        return (a, b)

    if isinstance(b, int):
        b = (b,)

    if len(a.shape) == 1:                # 1D array is a special case
        dim = b[-1]
        while a.shape[0] > dim:          # only shrink a
            if (dim - a.shape[0]) >= 2:  # trim off edges evenly
                a = a[1:-1]
            else:                        # or average adjacent cells
                a = 0.5*(a[1:] + a[:-1])
    else:
        for dim_idx in range(-(len(a.shape)),0):
            dim = b[dim_idx]
            a = a.swapaxes(0,dim_idx)        # put working dim first
            while a.shape[0] > dim:          # only shrink a
                if (a.shape[0] - dim) >= 2:  # trim off edges evenly
                    a = a[1:-1,:]
                if (a.shape[0] - dim) == 1:  # or average adjacent cells
                    a = 0.5*(a[1:,:] + a[:-1,:])
            a = a.swapaxes(0,dim_idx)        # swap working dim back

    return a

#start=datetime.datetime(2012,1,1,0,0)
#start = datetime.datetime.utcnow()
#tidx = netCDF4.date2index(start,tvar,select='nearest') # get nearest index to now
tidx = -1
timestr = netCDF4.num2date(tvar[tidx], tvar.units).strftime('%b %d, %Y %H:%M')

zlev = -1  # last layer is surface layer in ROMS
u = nc.variables['u'][tidx, zlev, j0:j1, i0:(i1-1)]
v = nc.variables['v'][tidx, zlev, j0:(j1-1), i0:i1]

lon_u = nc.variables['lon_u'][ j0:j1, i0:(i1-1)]
lon_v = nc.variables['lon_v'][ j0:(j1-1), i0:i1]
lat_u = nc.variables['lat_u'][ j0:j1, i0:(i1-1)]
lat_v = nc.variables['lat_v'][ j0:(j1-1), i0:i1]

lon=lon_rho[(j0+1):(j1-1), (i0+1):(i1-1)]
lat=lat_rho[(j0+1):(j1-1), (i0+1):(i1-1)]
mask = 1 - nc.variables['mask_rho'][(j0+1):(j1-1), (i0+1):(i1-1)]
ang = nc.variables['angle'][(j0+1):(j1-1), (i0+1):(i1-1)]

# average u,v to central rho points
u = shrink(u, mask.shape)
v = shrink(v, mask.shape)

# rotate grid_oriented u,v to east/west u,v
u, v = rot2d(u, v, ang)

from mpl_toolkits.basemap import Basemap
basemap = Basemap(projection='merc',llcrnrlat=39.5,urcrnrlat=46,llcrnrlon=-71.5,urcrnrlon=-62.5, lat_ts=30,resolution='i')
fig1 = plt.figure(figsize=(15,10))
ax = fig1.add_subplot(111)

basemap.drawcoastlines()
basemap.fillcontinents()
basemap.drawcountries()
basemap.drawstates()
x_rho, y_rho = basemap(lon,lat)

spd = np.sqrt(u*u + v*v)
h1 = basemap.pcolormesh(x_rho, y_rho, spd, vmin=0, vmax=1.0)
nsub=2
scale=0.03
basemap.quiver(x_rho[::nsub,::nsub],y_rho[::nsub,::nsub],u[::nsub,::nsub],v[::nsub,::nsub],scale=1.0/scale, zorder=1e35, width=0.002)
basemap.colorbar(h1,location='right',pad='5%')
title('COAWST Surface Current: ' + timestr);

Produces the following plot: enter image description here

[GIS] Plot in rstudio a NetCDF File – latbox by lonbox

1) It will be vastly easier if this works:

library(raster)
olr <- raster("t2m_50N70N.nc", varname = "var167")
tm2 <- brick("t2m_50N70N.nc", varname = "var179")

plot(olr[[1]])) ## first of nlayers(olr) bands

If this kind of mapping of the data is what you are after, I'd re-focus and work on learning raster. You'll also need ncdf4 package installed.

2) ncdf is no longer supported, you can now use ncdf4 as a safe default (raster does this for you).

3) To use NetCDF in the raw like this, you need to learn the start and count arguments to the "get var" function, in ncdf4 this is called ncvar_get. But, this is really painful since you have to juggle your time-subset across all the variables and work out from other variables what the right subset indexes for the three dimensions must be.

In a limited context, raster does everything for you - it's not general enough for all of NetCDF, but it's worth trying because if your files conform to its assumptions you get a lot for free.