I am working with a file which has the coordinates in WGS84 format. From this file I have calculated the Cartesian coordinates using the following code:

import numpy as np
import pyproj
from netCDF4 import Dataset
_projections = {}

def unproject(z, l, x, y):
    if z not in _projections:
        _projections[z] = pyproj.Proj(proj='utm', zone=z, ellps='WGS84')
    if l < 'N':
        y -= 10000000
    lng, lat = _projections[z](x, y, inverse=True)
    return (lng, lat)

def calculate_grid_definition(in_path):
    Calculate the latitude and longitude coordinates of the cell.

    :param in_path: Path to the file that contains all the information.
    :type in_path: str

    :return: Latitudes array, Longitudes array, Latitude interval, Longitude interval.
    :rtype: numpy.array, numpy.array, float, float
    with Dataset(in_path) as f:
        utm_lat = f.variables['LAT'][:]
        utm_lon = f.variables['LON'][:]
        geolat = np.zeros((len(utm_lon), len(utm_lat)))
        geolon = np.zeros((len(utm_lon), len(utm_lat)))
        for x, u_lon in enumerate(utm_lon):
            for y, u_lat in enumerate(utm_lat):
                lon_now, lat_now = unproject(32, 'N', u_lon, u_lat)
                geolat[x, y] = lat_now
                geolon[x, y] = lon_now

        return geolat, geolon

inpfile = 'path/to/filename.nc'
geolat, geolon = calculate_grid_definition(inpfile)

Now I also want to calculate the corner points of the grid cells. Assuming a regular rectangular grid, it would be very simple with the following code:

def create_bounds(coordinates):
    interval_x = coordinates[1, :] - coordinates[0, :]
    interval_y = coordinates[:, 1] - coordinates[:, 0]
    interval_x = np.tile(interval_x, (len(interval_y), 1))
    interval_y = np.tile(interval_y, (interval_x.shape[-1], 1)).T

    coords_left = coordinates - interval_x / 2
    coords_right = coordinates + interval_x / 2
    coords_bottom = coordinates - interval_y / 2
    coords_top = coordinates + interval_y / 2

    bound_coords = np.dstack((coords_left, coords_top,
                              coords_right, coords_bottom))
    return bound_coords

lat_bnds = create_bounds(geolat)
lon_bnds = create_bounds(geolon)

But, this is not right because the grids are skewed and result in an overlap. Any leads? For your information, you could download the original grid coordinates from this file. An alternative download link is here.

  • Access denied to file content :-( – wingnut May 14 at 0:49
  • Don't know why, but I have added an alternative download link from Onedrive – Vinod Kumar May 14 at 5:45

Are the original .nc coordinates in long/lat or grid coordinates? I ask, as it says LAT and LONG, rather than North and East. I guess they are UTM North and East coordinates in Zone 32 N, ie EPSG:32632.

If your original data are on a rectangular grid, then your unprojected data will not be.

Your function create_bounds makes boxes around the original UTM pixel centres after backprojection to geodetic coordinates.

Your interval_x etc computation uses specific rows of coordinates, but we know that the grid interval is not constant in geodetic coordinates. So you may be getting overlapping or disjoint polygons. Instead, maybe something like this would respect each individual grid element spacing:

dx = coordinates[1:,:] - coordinates[:-1,:]
dy = coordinates[:,1:] - coordinates[:,:-1]

Then the inter-row and -column differences would be retained. You could then compute the bounds of each cell from these.

Note that dx and dy outlined above will have M-1 X N-1 points, whereas you want M+1 X N+1 points to create bounds over all pixels. Just copy the edge values to one more row or column at each edge.

I looked at scipy.interpolate, but I couldn't see a quick way to generate midpoints and add edges.

  • Thanks a lot. You guessed it right, the names in the original files are wrong. These are north and east. Probably I could not convey what I want at the end. What you suggested provides the corner points f my entire data. What I actually want is corner points of the individual grids. For example, after transformation, geolat and geolon are array of dimension 642x873. The corner points should be array of shape 642x873x4. – Vinod Kumar May 13 at 16:40

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.