4

I am working on trying to get values for a specific latitude and longitude from a netcdf3 file with Python. The model uses X and Y coordinates and an LCC projection.

The netcdf file comes with a set of corresponding 2d lats and longs. I can look up my desired lat and long and find the nearest value of each in the 2d lat long file. I am using scipy.io netcdf and numpy. Dataset is a netcdf object.

lat = 41.050
lon = -125.685
wind_lats = dataset.variables['lat'][:]
wind_lons = dataset.variables['lon'][:]
lat_idx = np.abs(wind_lats - lat).argmin()
lon_idx = np.abs(wind_lons - lon).argmin()
lat_idx = np.unravel_index(lat_idx, wind_lats.shape)
lon_idx = np.unravel_index(lon_idx, wind_lons.shape)
wind_lat_y = lat_idx[0]
wind_lat_x = lat_idx[1]
wind_lon_y = lon_idx[0]
wind_lon_x = lon_idx[1]

This gives me a separate XY pair for the latitude and for the longitude. The variable of the netcdf file I am trying to access has indices time, height, y, and x.

It only uses 1 value for x and 1 for y. I have two values, but I don't know how to go from 2 to 1. Also the netcdf x and y variables are each a 1d array containing the respective x and y values.

The lat_idx and lon_idx after they are unraveled return the indices for the x and y values from the netcdf file. There is about a 100km difference between the x and y values of the lat longs. For example the latitude x value is -2641.2786. The longitude x value is -2543.75. The y values are similarly separated. The x and y values are not actual latitudes and longitudes themselves. The 2d lat long files reference the x and y values. Below is an example of how each of the files is organized:

From Panoply

2d latitude file:

y,          x,           lat
1959.0411,  -3043.5815, 38.9718
1959.0411,  -3031.3906, 38.9985
1959.0411,  -3019.1997, 39.0251

1d x file:

x,  x
-3043.5815, -3043.5815
-3031.3906, -3031.3906
-3019.1997, -3019.1997

1d y file:

y,  y
1959.0411,  1959.0411
1971.2321,  1971.2322
1983.4230,  1983.4231
  • How do you extract the two values from the netcdf file? It would be informative to know where they came from. How far apart are they? – Michael Stimson Nov 21 '17 at 1:58
  • Thanks for your time. I have updated the post with more info. – Anrydu Nov 21 '17 at 2:21
  • Those numbers don't look right.. values should be in the range -180 to 180 for X and -90 to 90 for Y. If you're accessing them by script can you post the lines of code that unravels these values, there might be a problem in your files or your script. – Michael Stimson Nov 21 '17 at 2:27
  • I am a bit confused. The x and y values are not themselves latitude and longitude values. The code to unravel the values is listed above. I find the nearest value to my desired latitude and longitude values from the netcdf lat and long arrays. Lat_idx and lon_idx unravel the indices returned from the previous argmin operations returning an x index and a y index. I end up with two sets of xy coordinates. One for lats and one for longs. – Anrydu Nov 21 '17 at 2:47
  • It seems that the values are not lat/lon degrees, but LCC-projected coordinates in meters or feet. You need to konw the parameters of the projection to convert them to real-world-lat/lon degrees. Usually the netcdf metadata contains all relevant information. – AndreJ Nov 21 '17 at 15:08
1

You need to loop over both dimensions x and y and compare the distance of every lon-lat pair in the x-y grids from your desired lon-lat point.

You seem to have indicated that in addition to the projection definition there are also supplemental 2D lon and lat arrays in the file, so you can just try to work with those rather than code up the map projection transformation equations.

So writing this very crudely in Java, it might look something like:

double ixmin = -1; // x index at closest grid point
double iymin = -1; // y index at closest grid point
double mindist = Double.INFINITY;

for (int iy = 0; iy < ylength; ++iy) {
    for (int ix = 0; ix < xlength; ++ix) {
        double dlon = mylon - lon(iy,ix);
        double dlat = mylat - lat(iy,ix);
        double dist = Math.sqrt (dlon * dlon + dlat * dlat);

        if (dist < mindist) {
            mindist = dist;
            iymin = iy;
            ixmin = ix;
        }
    }
}

System.out.println ("Closest grid point is at iy=" + iy + ", ix=" + ix);

Note that the above does not account for possible differences in handling the lon range. Some files will specify longitudes running from (0,360), others (-180,180). So you need to make sure that you don't have a dlon value whose abs value is over 360.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.