[GIS] Finding correct X Y values from two pairs of X Y map coordinates

coordinateslatitude longitudexy

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

Best Answer

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.

Related Question