Search code examples
pythonnetcdfpython-iris

how to convert projection x and y coordinate in netcdf iris cube to lat lon

see question above, just want to get lat/lon for Google maps from Iris cube...

cube.coord_system()

LambertAzimuthalEqualArea(latitude_of_projection_origin=54.9, longitude_of_projection_origin=-2.5, false_easting=0.0, false_northing=0.0, ellipsoid=GeogCS(semi_major_axis=6378137.0, semi_minor_axis=6356752.314140356))

enter image description here

Solution

  • If you only want to transform x- and y-coordinates to plot the data, you can use iris and cartopy:

    import iris
import numpy as np

    First, get coordinate points in the native projection

    proj_x = cube.coord("projection_x_coordinate").points
proj_y = cube.coord("projection_y_coordinate").points

    Next, make a pair of 2D arrays with the same shape

    xx, yy = np.meshgrid(proj_x, proj_y)

    Then extract the native projection and convert it to a cartopy projection:

    cs_nat = cube.coord_system()
cs_nat_cart = cs_nat.as_cartopy_projection()

    Next, create a target projection, e.g. the standard ellipsoid projection

    cs_tgt = iris.coord_systems.GeogCS(iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS)
# Again, convert it to a cartopy projection
cs_tgt_cart = cs_tgt.as_cartopy_projection()

    And finally, use cartopy's transform method to convert the 2D coordinate arrays in native projection to coordinates in the target projection.

    lons, lats, _ = cs_tgt_cart.transform_points(cs_nat_cart, xx, yy).T
# Note the transpose at the end.

    Also note that the function above always returns z coordinate array, but in this case is 0.

    You can then use lons and lats to plot your data in Google Maps or other applications. Remember that these new coordinates are curvilinear, so in practical terms have to be 2D arrays.

    If you, however, want to plot the data in iris (and matplotlib), you can do it this way:

    import cartopy.crs as ccrs
import matplotlib.pyplot as plt


proj_x = cube.coord("projection_x_coordinate").points
proj_y = cube.coord("projection_y_coordinate").points
cs_nat_cart = cube.coord_system().as_cartopy_projection()

fig = plt.figure()
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.pcolormesh(proj_x, proj_y, cube.data, transform=cs_nat_cart)
ax.coastlines()

    Hope this helps.