examples outdated?

[boristyukin]

hey guys, really neat architecture and documentation! very impressive work!

I noticed that MOGREPS-UK files now have projected x/y coordinates vs. lat/lons. Being new to all of this, I have hard time using projected x/y. I was only able to figure out that recent files are using iris.coord_systems.LambertAzimuthalEqualArea

can you tell me how to get regular lat/lons from MOGREPS-UK so I can use them with Google Map? I need to extract data from cubes and store it in another system in a simple table for my needs.

For example, one of projected y/x in a cube is -158000.0,-36000.0 but google maps does not recognize this as valid coordinate.

I understand I need to change projection of the cube but have no clue now.

Would be great to update examples as well.

[AlexHilson]

Hi, thanks for the comments. I'm afraid I don't have a perfect example to hand for how to do this, but can offer some links. If you have success in this and can share code we'd love to add examples for doing this reprojection, I think it's a question a lot of people will have.

It's good that Iris has picked up the co-ordinate system (reference https://scitools.org.uk/iris/docs/latest/iris/iris/coord_systems.html#iris.coord_systems.LambertAzimuthalEqualArea). That class has a 'as_cartopy_crs' method. Cartopy is commonly used in the Iris ecosystem for handling this kinds of transformations.

The Cartopy CRS reference is found at https://scitools.org.uk/cartopy/docs/latest/crs/index.html#cartopy.crs.CRS. So what might work for you is to create a Mercator CRS (https://scitools.org.uk/cartopy/docs/latest/crs/projections.html#mercator), then use the transform_points method of the CRS class to convert from LambertAzimuthalEqualArea to Mercator. That would give you a lat/lon value for that point (I'm not sure what projection GoogleMaps is looking for, I suspect Mercator will be ok but you may need to translate again into WebMercator).

One thing to be wary of is that the correct properties for LambertAzimuthalEqualArea have been picked up from the source file, i.e. that properties like latitude_of_projection_origin and longitude_of_projection_origin are set rather than using the default of 0. I believe the correct values for the center would be longitude -2.5, latitude 54.9 but I hope that the correct values are included in the file. (edit: I believe false eastings / northings should be 0, and the ellipsoid is GRS80)

Just for reference in case Cartopy doesn't work for you, Cartopy is using proj under the hood. Any set of proj tooling or libraries would be able to translate between lambert azimuthal co-ordinates and lat/lon co-ordinates (given the correct definition of both reference systems). E.g. I think the invproj command line tool from https://proj.org/apps/proj.html would do what you're looking for.

[AlexHilson]

Looked around and I believe that this is the kind of command you want

invproj +proj=laea +lon_0=-2.5 +lat_0=54.9 +x_0=0.0 +y_0=0.0 +ellps=GRS80 <<EOF
-36000 -158000
-1158000 902000
EOF
3d2'31.982"W    53d28'45.337"N
24d30'35.729"W  61d19'7.892"N

[boristyukin]

@AlexHilson thanks for the hints! I will play a bit more. I never worked with this sort of data so it is a bit challenging.

this what cube.coord_system() return

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))

[boristyukin]

BTW I tried this iris cube projection described here https://scitools.org.uk/iris/docs/v1.8.0/whatsnew/1.0.html#cube-projection

but it seems it does not like MOGREPS-UK coordinates and I see some hardcoded lines expecting lat/lon coordinates specifically

~/miniconda3/envs/venv/lib/python3.8/site-packages/iris/analysis/cartography.py in project(cube, target_proj, nx, ny)
    595         lon_coord, lat_coord = _get_lon_lat_coords(cube)
    596     except IndexError:
--> 597         raise ValueError('Cannot get latitude/longitude '
    598                          'coordinates from cube {!r}.'.format(cube.name()))
    599 

ValueError: Cannot get latitude/longitude coordinates from cube 'air_temperature'.

[boristyukin]

tried your idea with cartopy but it did not work. produced some crazy numbers, certainly not degrees or anything like that.

I am pulling my hair out. I need to solve it in python, cannot use command line tool for that (it is part of larger code)

import cartopy.crs as ccrs
target_proj = ccrs.Mercator(central_longitude=0, false_easting=0.0, false_northing=0.0)
source_proj = cube.coord_system().as_cartopy_crs()
source_manual = ccrs.LambertAzimuthalEqualArea(central_longitude=-2.5, central_latitude=54.9, false_easting=0.0, false_northing=0.0, globe=cartopy.crs.Globe(ellipse='GRS80'))
print(target_proj) # target crs
print(source_proj) # source crs
print(target_proj.transform_point(-158000.0,-36000.0, source_proj))
print(target_proj.transform_point(-158000.0,-36000.0, source_manual))

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)) <cartopy.crs.Mercator object at 0x7fa50e6f1540> <cartopy.crs.LambertAzimuthalEqualArea object at 0x7fa5196ad180> (-550185.3208803544, 7240551.325457196) (-550185.3208803543, 7240551.325457192)

[boristyukin]

I gave up on cartopy, it just does not make sense.

But this below worked for me! thanks so much @AlexHilson for the idea and command-line syntax. I would take me days to figure it out!

import pyproj
p = pyproj.Proj("+proj=laea +lon_0=-2.5 +lat_0=54.9 +x_0=0.0 +y_0=0.0 +ellps=GRS80")
# lon, lat = p(x, y, inverse=True)
lon, lat = p(-36000, -158000, inverse=True)
print (lat, lon)

53.47926025120268 -3.0422171310489534

[AlexHilson]

Thanks @boristyukin. I'm going to open a new issue referencing the code you show above so we move towards having a full example for working with LAEA projected data. Sorry to hear Cartopy didn't work for you, if you find yourself looking at it again you might have luck asking on Stack Overflow with Iris / Cartopy tags.