iris.coord_system for EPSG:3857 or Pseudo-Mercator

[ecasellas]

Hi,

I'm using iris together with IMPROVER to obtain a blending of probability of precipitation forecasts. Today I was struggling with a regridding problem, which I think it is related with the definition of the coordinate system of the cube that defines the target grid.

The thing is that I have a Cube which coordinates have iris.coord_system.TransverseMercator() as coordinate system and I would like to regrid it to a Cube with EPSG:3857 or Pseudo-Mercator projection. However, I'm failing to obtain a proper coord_system for this projection and I don't know if it's supported by any of the coordinate systems available in iris (maybe a special case of iris.coord_system.Mercator()).

To solve this I was thinking of transforming the iris cube to xarray, then using rio.reproject('epsg:3857'), and finally processing the output. Still, I was wondering if there is any way to do this with iris only; maybe the solution is very simple and I'm wrongly defining some parameters, but I don't get it.

Thank you very much for your time,

Enric

[pp-mo]

Thanks @ecasellas for engaging !

I haven't just "solved" this, but I have some suggestions to maybe cast some light on it ...

I looked here for info.
It says there, on the one hand, "Remarks: Not a recognised geodetic system. Uses spherical development of ellipsoidal coordinates. Relative to WGS 84 / World Mercator (CRS code 3395) gives errors of 0.7 percent in scale and differences in northing of up to 43km in the map (21km on the ground)."
But on the other hand, it gives proj prarams of :
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
So, that looks pretty standard to me, and it seems to imply that it is basically a plain Mercator projection (not transverse or anything odder).

So on the face of it, it looks like it might simply be the same as a plain cartopy.crs.Mercator() or iris.coord_systems.Mercator(), with the appropriate spherical datum,
So I'd say it is certainly worth trying those, if you have suitable data, to see if a data plot lines up with coastlines as you would expect.

But I do have remaining doubts, as there are things here which I still don't understand ..

• firstly from the "remarks;" comment above, what is meant by "spherical development of ellipsoidal coordinates" ?
• and, in the above proj4 params, what the +nadgrids=@null part means,
  • I couldn't find any specific info on that "nadgrids" parameter, so we are in realms of speculation.
    -- though there is info here that might explain it, if properly understood (but I failed to do so!)
  • however, I have a thought that it might be specifying a specific ellipsoid/datum (which could be consistent with the 'remarks')
  • also, it should be fairly simple to at least visually compare plots based on cartopy crs with + without this argument : Since we can now create a Cartopy projection directly from specific PROJ parameters.
    . * ... or for that matter, to compare coordinate values transformed from another system, like plain lat-lon or your

Possibly... what it all means is that the projection combines "real" geographical lats+lons with a spherical datum, leading to certain distortions.

Can anyone cast any more light on those things ?

[ecasellas]

Hi @pp-mo,

Thank you very much for your rapid response! I tried to define an iris.coord_systems.Mercator() with the proj params you indicated (I could not inform all of them), but still, the reprojection doesn't look good. I have to admit that this is a bit of a headache for me, but I think that I found something that may explain why I can't properly define Pseudo-Mercator projection with iris.coord_systems.Mercator().

Searching in cartopy source code, I found the definition of Mercator.GOOGLE which is an equivalent name to Pseudo-Mercator or Web-Mercator, here. ccrs.Mercator in cartopy allows defining some parameters that, I believe, are not available in iris.coord_systems.Mercator(): min_latitude, max_latitude, and part of globe.

Mercator.GOOGLE is defined in cartopy as

# Define a specific instance of a Mercator projection, the Google mercator.
Mercator.GOOGLE = Mercator(min_latitude=-85.0511287798066,
                           max_latitude=85.0511287798066,
                           globe=Globe(ellipse=None,
                                       semimajor_axis=WGS84_SEMIMAJOR_AXIS,
                                       semiminor_axis=WGS84_SEMIMAJOR_AXIS,
                                       nadgrids='@null'))

Mercator in iris is

class Mercator(CoordSystem):
    """
    A coordinate system in the Mercator projection.

    """

    grid_mapping_name = "mercator"

    def __init__(
        self,
        longitude_of_projection_origin=None,
        ellipsoid=None,
        standard_parallel=None,
        scale_factor_at_projection_origin=None,
        false_easting=None,
        false_northing=None,
    ):

And in iris the _globe method is

    @cached_property
    def _globe(self):
        """
        A representation of this CRS as a Cartopy Globe.

        Note
        ----
        This property is created when required and then cached for speed. That
        cached value is cleared when an assignment is made to a property of the
        class that invalidates the cache.
        """
        if self._datum is not None:
            short_datum = _short_datum_names.get(self._datum, self._datum)
            # Cartopy doesn't actually enact datums unless they're provided without
            # ellipsoid axes, so only provide the datum
            return ccrs.Globe(short_datum, ellipse=None)
        else:
            return ccrs.Globe(
                ellipse=None,
                semimajor_axis=self._semi_major_axis,
                semiminor_axis=self._semi_minor_axis,
            )

So I believe that (correct me if I'm wrong), at the moment, it's not possible to define nagrids as '@null', may this be critical for the projection. I also admit that I don't fully understand how it works.

I think I haven't answered any of the remaining doubts you wrote, but maybe we have now more information in the discussion regarding the possibility to use Pseudo-Mercator in iris compared to cartopy. Indeed, if EPSG:3857 is not recognised as a geodetic system, it may not be included in iris. I raised the question to try to only use iris in the software we are developing, basically in the part of obtaining data in EPSG:3857 for visualisation purposes in a Leaflet map. There are some workarounds to obtain these data, so I don't want this to become a waste of time for you :). However, I found very interesting the questions you raised, and it would be nice if someone cast more light on them.

Again, thank you very much!

[pp-mo]

Now you point out what is in Cartopy, it seems clear to me that the "+nadgrid=@null" is the key thing that make a difference
-- and I guess it probably does mean something like "no datum - use spherical earth" (i.e. it sets something which determines the datum to 'null' i.e. nothing ??).

I'm frankly a bit surprised that this makes much visible difference to a plot, but you are saying it creates visible problems ??

[pp-mo]

More to the point ... I think I see more clearly now what you are asking for.

The point is, that the range Iris coordinate systems is simply more limited than Cartopy, because we it supports what can be expressed in the netcdf-CF encoding, so that is limited to these ones. And as you can see there, their support for different earth ellipsoids is rather limited.

Of course you don't need to use Iris for plotting : You can use Cartopy directly, which now gives you access to any PROJ projection.
So, you should be able to make plots + check if it is making correct sense of the data.
This may help : https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html?highlight=project
Iris is using this internally -- the Iris plotting routines are only a fairly thin wrapper to Cartopy (mostly, they just work out which coordinates to plot with).

Regridding is a rather different problem, though.
I guess you ideally would want to use Iris for that, but I don't think we will have anything which obviously fits the requirement.
( Unfortunately I know nothing about regridding in xarray, so I can't say whether that is likely to work effectively ).

A really key question is, what kind of regridding ?
For conservative regridding (i.e. weighting by areas overlap between input and output gridcells), I'm fairly confident that we don't have a solution which can work with cells defined in an arbitrary coordinate system, since the geometric calculations are very involved.
For nearest-neighbour regridding, I think the calculation certainly is possible, and I think maybe the UnstructuredNearestNeighbour scheme can do it, but it is probably very inefficient -- and the data needs to be transformed in some rather non-obvious ways to make it work.

[pp-mo]

Call out to @stephenworsley ...
are you aware of any approach in Iris or iris-esmf-regrid that might be able to help with this?

The target will be a structured "Cube with EPSG:3857 or Pseudo-Mercator projection", and the source is a structured cube with plain 'TransverseMercator' grid -- but I'm thinking that the source can be transformed into an unstructured form, with the TM gridpoints projected into the target coordinate system (i.e. a different pseudo-mercator X and Y at every point).

I'm fairly convinced that you can do this calculation also by presenting the data to UnstructuredNearestRegridder, since it is based on a trajectory regridding approach.
But it's a fuss to re-form the input in line with the requirements of that, and it doesn't provide weights caching.
Basically there is a potential solution in there trying to get out, based on the KDtree trajectory regridding, but it is not made fully available in the user API.

[ecasellas]

Hi to all,

First of all, I think I found a problem in my code (I'm feeling a bit ashamed right now). The source cube did not have a proper ellipsoid definition, because I used semi_major_axis and inverse_flattening parameters, but I didn't use the named arguments properly.

ellipsoid=iris.coord_systems.GeogCS(proj['semi_major_axis'],
                                                proj['inverse_flattening']))

So the coord_system assumed semi_minor_axis as inverse_flattening. Now, defining source ellipsoid as

ellipsoid=iris.coord_systems.GeogCS(semi_major_axis=proj['semi_major_axis'],
                                                inverse_flattening=proj['inverse_flattening']))

it worked well.

In addition, I modified Iris source code to include in coord_systems.Mercator() the arguments min_latitude and max_latitude, and also modified GeogCS to include the parameter nadgrids='@null' when defining the equivalent cartopy Globe. The thing is that this produced no change in the reprojection process, and it is, as far as I could see, the same result. Still, I don't know if this should be the correct behaviour accounting to the differences between Mercator and Pseudo-Mercator, but in the end, it worked. I will try to investigate a bit further why the behaviour is the same in both cases.

Thank you very much for all your answers and options you offered me, it also make me discover some new iris features such as the UnstructuredNearestNeighbour or iris-esmf-regrid.

[stephenworsley]

@pp-mo In terms of being able to understand different projections, iris-esmf-regrid should be able to recognise the same kinds of projections that iris can, since it is built off iris objects. If iris is has a problem with the EPSG:3857 projection then i expect iris-esmf-regrid might also. With that said, if you find a way to transform your projected coordinates into 2D lat-lon coordinates iris-esmf-regrid might be able to accomplish the rest of the regridding task. At the moment, I believe only area weighted conservative regridding is available via iris-esmf-regrid (though bilinear and unstructured nearest should hopefully be added soon).

[pp-mo]

Thanks for explaining @ecasellas 💐
I just wrote a reply, then deleted on re-reading yours, which I think answered all my questions !

So... are you actually ok with this now, should we close the issue ?

[ecasellas]

Yes, sure! Thank you very much @stephenworsley and @pp-mo!