Plots of data with a CRS without EPSG code not possible

[timschmi95]

Calculating impacts with a CRS not corresponding to an EPSG code (e.g. a rotated lat/lon coordinate grid with a specified center point: ccrs.RotatedPole(pole_latitude=30, pole_longitude=-175)) works well, but none of the plotting functions (hazard, exposure, impact) work.
The issue is that u_plot.get_transformation() is called, which tries and fails to get an EPSG code from the CRS.
For most plots, simply directly passing the CRS directly instead of trying to get an EPSG code in the u_plot.get_transformation() function creates a correct plot.

Minimal Rep. Example for (failing) exposure plots:

import geopandas as gpd
from shapely.geometry import Point
from climada.entity import Exposures

#Define a rotated pole projection
prj = ccrs.RotatedPole(pole_latitude=30, pole_longitude=-175)

# Create grid points (2x4 grid for 8 points)
x_coords = [0, 1, 2, 3]
y_coords = [-2,-1,0, 1]
points_list = [Point(x, y) for y in y_coords for x in x_coords]

# Random values
values = np.random.rand(16)

# Create GeoDataFrame
gdf = gpd.GeoDataFrame({'value': values}, geometry=points_list, crs=prj)

# Create Exposures object and plot data
exp = Exposures(gdf)
ax=exp.plot_raster()

This example results in a rasterio.errors.CRSError from trying epsg = CRS.from_user_input(crs_in).to_epsg() in the get_transformation() function.

My suggestion would be to catch the CRSError if no EPSG code can be found, and in this case directly pass the input crs (while checking it's the correct object type and potentially other checks).
Is this suitable, or is there something speaking against allowing plots for data in a CRS without an EPSG code?

[peanutfun]

The entire get_transformation seems superfluous to me. But this might be a recent development. The GeoDataFrame.crs is a pyproj.CRS instance, and one should be able to create a cartopy.crs.CRS from that via

gdf = geopandas.GeoDataFrame(..., crs=...)
projection = cartopy.crs.CRS(gdf.crs)

See SciTools/cartopy#813 (comment)

Your code is a bit odd, though. Cartopy CRS are not directly compatible to PyProj CRS and intended for plotting only, at least in my understanding. I don't understand why you use it to instantiate a GeoDataFrame. Note also that cartopy distinguishes between the projection of the plot and the projection of the data. In general, you should do something like this:

# Here, you set the projection your data is in.
# You cannot choose it, but you can reproject, of course
gdf = geopandas.GeoDataFrame(..., crs="EPSG:4326")

# Here, you choose the projection of your plot
fig, ax = plt.subplots(subplot_kw={"projection": cartopy.crs.Mercator()})

# Here, you tell cartopy what projection the data is in
# (if you don't it assumes the same as the plot)
ax.plot(gdf[column], transform=cartopy.crs.CRS(gdf.crs))

[timschmi95]

I agree, it would be neat not to always try to find the EPSG, but directly use the CRS from the GeoDataframe (respectively the Exposures.crs attribute). If this is always a pyproj.CRS (or cartopy.crs) instance, that shouldn't be a problem. I didn't suggest that as I wasn't sure if it could also be e.g. a proj4 string or wkt string, which is currently handled with the pyproj.CRS.from_user_input() method in the get_transformation() function.

Regarding the use of cartopy.crs as GeoDataFrame.crs: My data are defined on a RotatedPole coordinate system, and I only found this in the cartopy.crs list, not in the pyproj.CRS. Of course, it can be reprojected, but then wouldn't be on a regular grid anymore.

(In general; outside climada) I commonly use cartopy.crs (which can be geographical cartopy.crs.CRS or projected cartopy.crs.Projection) for both the data coordinates and the map projection coordinates. This usually works well for me, whether I use the same map projection as the data or a different projection. Do you think this is discouraged? I.e. would you not want to allow the use of GeoDataframes that have a Cartopy CRS instead of a PyProj CRS as the GeoDataFrame.crs attribute?

[peanutfun]

Interesting! I haven't seen the use of cartopy CRS in other circumstances than plotting.

If this is always a pyproj.CRS (or cartopy.crs) instance, that shouldn't be a problem.

It is according to the geopandas docs 🤔

Does this code work in your case? (specifically, the cartopy.crs.CRS(gdf.crs)) part)

# Here, you choose the projection of your plot
fig, ax = plt.subplots(subplot_kw={"projection": cartopy.crs.Mercator()})

# Here, you tell cartopy what projection the data is in
# (if you don't it assumes the same as the plot)
ax.plot(gdf[column], transform=cartopy.crs.CRS(gdf.crs))

[chahank]

Interesting! I haven't seen the use of cartopy CRS in other circumstances than plotting.

It is used for re-projection or area calculations for instance.

I.e. would you not want to allow the use of GeoDataframes that have a Cartopy CRS instead of a PyProj CRS as the GeoDataFrame.crs attribute?

Would the best maybe be to let geopandas handle the crs?

[timschmi95]

Does this code work in your case? (specifically, the cartopy.crs.CRS(gdf.crs)) part)

fig, ax = plt.subplots(subplot_kw={"projection": cartopy.crs.Mercator()})

ax.plot(gdf[column], transform=cartopy.crs.CRS(gdf.crs))

The function ax.plot does not lead to a geographical plot for me (but neither if I initialize the dataframe with "EPSG:XXXX" and the .crs attribute is a pyproj.crs.CRS instance. But the equivalent with a scatter plot works also well with a cartopy.crs as the crs of the geodataframe:

fig, ax = plt.subplots(subplot_kw={"projection": ccrs.Mercator()})
p1= [ccrs.RotatedPole(pole_latitude=30, pole_longitude=-175),[-2,0,2,4],[-2,0,2,4]]
prj, x_coords, y_coords = p1 # or p2 or p3
points_list = [Point(x, y) for y in y_coords for x in x_coords]

gdf = gpd.GeoDataFrame({'value': values,'impf_HL':1}, geometry=points_list, crs=prj)
ax.scatter(gdf.geometry.x, gdf.geometry.y, transform=ccrs.CRS(gdf.crs), c=gdf['value'])

[peanutfun]

Great! This means we can rely on the ccrs.CRS(gdf.crs) instead of calling get_transformation.

I also noticed that all Climada code actually ignores the units returned by get_transformation, so it seems safe to remove the function altogehter (see https://github.com/search?q=repo%3ACLIMADA-project%2Fclimada_python%20get_transformation&type=code)

[timschmi95]

One point to consider is that the get_transformation() function also converts EPSG:4326 to PlateCarree() and EPSG:3395 to Mercator().

Current code in get_transformation():

    try:
        epsg = CRS.from_user_input(crs_in).to_epsg()
        if epsg == 3395:
            crs = ccrs.Mercator()
        elif epsg == 4326:  # WSG 84
            crs = ccrs.PlateCarree()
        else:
            crs = ccrs.epsg(epsg)
    .......

And the returned crs is then not just used as the data coordinates, but also as the projection coordinates for the map (e.g. in the u_plot.make_map() function). Which is why the conversion is necessary.

So if we remove the function, we would need to do some edits in the downstream plotting function to make sure that the ccrs.CRS(gdf.crs) is only used as map projection if it is in fact possible (i.e. add the same two exceptions that convert the 2 aforementioned EPSG codes to their corresponding map projection). I think most other EPSG-code-CRS can directly be used as projections as well (at least all the ones I've encountered, and apparently also all the ones other climada users have used, since no issues were raised regarding it)

Alternatively (potentially cleaner), we can keep the function, but remove the return of units and simply pass ccrs.CRS(crs_in) if it is not EPSG:4326 or EPSG:3395. This would also allow keeping projections that are not EPSG codes.

[peanutfun]

I'm not sure I understand. PlateCarree (as transform parameter) is the cartopy equivalent to WGS84, so I would hope that
ccrs.CRS(proj.crs.CRS.from_user_input("EPSG:4326")) and ccrs.PlateCarree() result in the same thing. Same goes for EPSG:3395, which is a global Mercator projection.

[timschmi95]

The difference is quite subtle, but it is not possible to initialize a map e.g. with

proj = ccrs.CRS(CRS.from_user_input("EPSG:4326")) 
fig,ax = plt.subplots(1,1,subplot_kw={'projection':proj})

one instead needs to use PlateCaree().
However, the coordinates are equivalent (both with WGS84). But only PlateCarree() is a 'projected' coordinate system, the EPSG:4326 lacks the definition of how exactly to project the (same) coordinates onto a 2D map.

Also note that apart from plotting (e.g. for simply transforming points from some third CRS, both PlateCarree() and ccrs.CRS(proj.crs.CRS.from_user_input("EPSG:4326")) work.

[peanutfun]

That's correct! I think the subtlety here is that what I propose only works in ax.plot(..., transform=crs), whereas one really has to use a cartopy CRS for defining the plot projection (as in your example). However, get_transformation() is only used for the former. Or do I see this wrong? 🤔

[timschmi95]

The returned value of get.transformation() is used for both, not only as data transform. It ends up as transform argument e.g. in ax.pcolormesh (for the data coordinates), but it is also used in the make_map function (i.e. to define the projection).
E.g. in the Hazard.plot_intensity() function, the returned value of get.transformation() is eventually passed to u_plot.geo_im_from_array where it is used for the map projection as well.

[peanutfun]

Thanks for digging that up! 🙌 It might make sense to separate the two. I'll have a look and come up with a proposal