How to properly set projection and transformation in cartopy geoaxes in matplotlib plotting

Question

I am facing serious difficulties in using Python geopandas, cartopy and matplotlib to work together in a proper plot of my shapefile data.

The issue comes from the difficulty in properly setting the Transform and the Projection objects of my shapefile data.

The example I am describing here is relative to a SHP in SIRGAS 2000 projection, whose WKT format is:

GEOGCS["SIRGAS 2000",DATUM["D_SIRGAS_2000",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

From the WKT above, one may see that my coordinates are in degrees. The globe is not the WGS84, but GRS_1980, which is similar (but not equal) to WGS84.

This crs is not yet implemented in cartopy. So trying to use functions as cartopy.crs.EPSG(4674) doesn't work, especially because this is not a projection crs, but a geographic crs.

Following to the plotting Issue, I logically assumed that I should use one of the four cartopy.crs options to plot my data:

transform_option_1 = cartopy.crs.PlateCarree() # works well

transform_option_2 = cartopy.crs.PlateCarree(globe=ccrs.Globe(ellipse='GRS80')) # works well

transform_option_3 = cartopy.crs.Geodetic(globe=ccrs.Globe(ellipse='GRS80')) # doesn't work

# using wkt:

WKT = """GEOGCS["SIRGAS 2000",DATUM["D_SIRGAS_2000",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]"""

transform_option_4 = ccrs.PlateCarree(WKT) # results in an error message (see message in annex)

However, when I try the four transform options to plot my data, an error appears for my transform_option_3 and transform_option_4 options.

The first two options (transform_option_1, transform_option_2), everything works well. Also when I compare the first two Transform options in the plot, they result in the same figure, though they should result in slightly different plots (due to the globe object). So here is my first question:

Question 1: why the ccrs.PlateCarree() object result in the same plot for my data, despite different Globe objects?

Question 2: since cartopy accepts WKT string format for transform object instanciation, why an error appear? See error message annexed below this message (error message number 1).

Question 3: what is the difference between geodetic and platecarree projections in cartopy? I can't figure out when I should use or the other. Since the description of both projections imply coordinates in degrees, their difference is strange to me. Also, when I try both options for plotting my SHP using cartopy, only the PlateCarree works. The other results in a black geoaxes.

Here is a code snippet for better clarification of my Question 3:

# importing libraries:
import geopandas as gpd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs


# importing data:

SHP_path = 'Data_SHP.shp'

SHP = gpd.read_file(SHP_path)


# setting projection to PlateCarre:

projection = ccrs.PlateCarree() # projection.proj4_init


# setting SHP data transform

Transform2 = ccrs.Geodetic(globe=ccrs.Globe(ellipse='GRS80'))


# creating figure:

fig, ax = plt.subplots(1, subplot_kw={'projection':projection})


SHP.plot(ax=ax, transform=Transform2, 
         color='k', alpha=0.5)

# defining a small function to add grilines to the plot:

def custom_plot(geoaxes, projection):



    gl = geoaxes.gridlines(crs=projection)

    tick_axis_positions={'xlabels_top':False,
                     'ylabels_left':True,
                     'ylabels_right':False,
                     'xlabels_bottom':True}

    gl.xlabels_top = tick_axis_positions['xlabels_top']
    gl.ylabels_left = tick_axis_positions['ylabels_left']
    gl.ylabels_right= tick_axis_positions['ylabels_right']
    gl.xlabels_bottom = tick_axis_positions['xlabels_bottom']


    return geoaxes


# adding gridlines to the geoaxes:

ax = custom_plot(ax, projection)

fig.show()

# end of code

Question 4: Why can't I instanciate a cartopy Projection object using my crs (in WKT format)? When I try that, an error message appears (see error message 2 in annex).

Projection = ccrs.Projection(proj4_params=WKT, globe=ccrs.Globe(ellipse='GRS80'))

References:

1) discussion relating how to create cartopy.crs transform object from WKT or proj4 coordinate system projection formats

2) Mapping proj4 to cartopy CRS

3) Plotting a straight line in Cartopy, Robinson projection

4) example of plotting SHP with cartopy

5) cartopy transform and projection theory

---

Annex:

Error message 1)

TypeError: unsupported operand type(s) for -: 'str' and 'float'


Traceback (most recent call last):
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\backends\backend_qt5.py", line 519, in _draw_idle
    self.draw()
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\backends\backend_agg.py", line 402, in draw
    self.figure.draw(self.renderer)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\figure.py", line 1649, in draw
    renderer, self, artists, self.suppressComposite)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\mpl\geoaxes.py", line 385, in draw
    inframe=inframe)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\axes\_base.py", line 2628, in draw
    mimage._draw_list_compositing_images(renderer, self, artists)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\collections.py", line 262, in draw
    transform, transOffset, offsets, paths = self._prepare_points()
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\collections.py", line 240, in _prepare_points
    for path in paths]
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\collections.py", line 240, in <listcomp>
    for path in paths]
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\transforms.py", line 2451, in transform_path_non_affine
    return self._a.transform_path_non_affine(path)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\mpl\geoaxes.py", line 173, in transform_path_non_affine
    src_path.vertices, self.source_projection)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\crs.py", line 746, in quick_vertices_transform
    bboxes, proj_offset = self._bbox_and_offset(src_crs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\crs.py", line 709, in _bbox_and_offset
    lon_0_offset = other_lon_0 - self_lon_0
TypeError: unsupported operand type(s) for -: 'str' and 'float'
Traceback (most recent call last):
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\backends\backend_qt5.py", line 519, in _draw_idle
    self.draw()
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\backends\backend_agg.py", line 402, in draw
    self.figure.draw(self.renderer)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\figure.py", line 1649, in draw
    renderer, self, artists, self.suppressComposite)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\mpl\geoaxes.py", line 385, in draw
    inframe=inframe)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\axes\_base.py", line 2628, in draw
    mimage._draw_list_compositing_images(renderer, self, artists)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\image.py", line 138, in _draw_list_compositing_images
    a.draw(renderer)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\collections.py", line 262, in draw
    transform, transOffset, offsets, paths = self._prepare_points()
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\collections.py", line 240, in _prepare_points
    for path in paths]
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\collections.py", line 240, in <listcomp>
    for path in paths]
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\matplotlib\transforms.py", line 2451, in transform_path_non_affine
    return self._a.transform_path_non_affine(path)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\mpl\geoaxes.py", line 173, in transform_path_non_affine
    src_path.vertices, self.source_projection)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\crs.py", line 746, in quick_vertices_transform
    bboxes, proj_offset = self._bbox_and_offset(src_crs)
  File "C:\Users\lealp\AppData\Local\conda\conda\envs\Python_3.7\lib\site-packages\cartopy\crs.py", line 709, in _bbox_and_offset
    lon_0_offset = other_lon_0 - self_lon_0
TypeError: unsupported operand type(s) for -: 'str' and 'float'

Error message 2)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-165-9916b2ddcc67> in <module>
----> 1 ccrs.Projection(proj4_params=Transform_from_wkt, globe=ccrs.Globe(ellipse='GRS80'))

TypeError: Can't instantiate abstract class Projection with abstract methods boundary, threshold, x_limits, y_limits

Answer 1

1. PlateCarree is faking using lat/lons "projection" coordinates. What it is really is Equidistant Cylindrical with a special default globe set that allows lat/lons to work. See more in this bug. Because of this, I don't think passing a custom GRS80 globe to PlateCarree is actually working for you. That bug outlines it breaking everything.

2. I'm not aware of CartoPy actually accepting WKT anywhere. Creating a projection from a Proj.4 string or WKT I think is an open issue.

3. PlateCarree works by taking lat/lon as coordinates on a plane, so lines between lat/lon points are straight lines on a plane. It's a gross simplification, but it makes many things work, like contouring. Geodetic treats lat/lons properly as coordinates on a sphere, so lines connecting lat/lons are great circles. It's the proper way to do things, but it makes the math more complicated and won't work for some plotting options. See this example to see how it changes.

4. See (2) above

I'm curious what problems you ran into passing a custom ellipsoid to Geodetic. I had no problems with this code and Cartopy 0.17:

import cartopy.crs as ccrs
globe = ccrs.Globe(ellipse='GRS80')
crs = ccrs.Geodetic(globe=globe)