Plotting Spatial Vector Data in Projected Axes by Matplotlib (Python)

Dear all,

I am a remote sensing PhD student from the National Spatial Research Institute of Brazil. I am quite new to Python, specially using cartopy and geopandas for Map creations. Up till now, I have been using softwares as Quantum-GIS and ARC-GIS for my spatial analysis.

Since my migration to Python, I am facing serious difficulties in order to understand how one may plot a spatial vector data (i.e. Shapefiles) in to a geographic projected axes (geoAxes) of Matplotlib library from Python.

My question arrises from the fact that one may have to handle multiple projections in this attempt: one projection from the Shapefile itself; one projection for the creation of the Geoaxes; and lastly, one projection for inserting the SHP geometric data into the Geoaxes.

Here is the case scenario:

1) Assume a case Shapefile (SHP) opened as a GeoDataFrame by geopandas in Python. It has a certain projection "A" (i.e.: SIRGAS 2000, EPSG: 4674), for example.

2) The SHP projection (proj4 string or its EPSG code) is not implemented in Cartopy Projection list.

With that in mind, here is a code snippet that illustrastes my question:

Coding.....

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

    SHP = gpd.read_file('my_file.shp') 

    # SHP.crs.EPSG == 4674    --->> True

    CRS_Axes = ccrs.AlbersEqualArea()   # Projection 2. Note that it is different from the SHP crs!

    Fig, Axes = plt.subplots(nrows=1, ncols=1, subplot_kw={'projection': CRS_Axes})

    SHP.plot(ax=Axes, transform=ccrs.Geodetic()) # Projection 3. Also different from the other two CRSs

End of Code...

So here are my questions:

1) Why set to Geodetic projection (or any other projection that is not necessarily the SHP, nor the Axes CRS) in order to insert the SHP geometric data into the Geoaxes?

2) How does Matplotlib deals with that combination of possibly three different CRSs?

3) Does Matplotlib reproject the geometric data from SHP to the Geodetic system, and later on it reprojects into the Axes CRS? If so, why not reproject the geometry data from SHP directly into the Axes CRS? In this way, the programmer/user would not have to insert a third party projection into the code. Code example below:

Code example:

    CRS_Axes = ccrs.AlbersEqualArea()   # Projection 2. Note that it is different from the SHP crs!

    Fig, Axes = plt.subplots(nrows=1, ncols=1, subplot_kw={'projection': CRS_Axes})

    SHP.plot(ax=Axes, transform=CRS_Axes) # Projection 2. The same projection from the Geoaxes

End of Code

4) Last possible scenario would be reproject the SHP into the GeoAxes CRs prior to the plotting. This is the only possible case, in which I believe one can be certain of the reprojected coordinates of the SHP geometric data and the GeoAxes grid. Code example below:

Code example

       CRS_Axes = ccrs.AlbersEqualArea()

       CRS_Axes_proj4 = CRS_Axes.proj4_init

       # reprojecting SHP into the same CRS that will be used in the GeoAxes creation:

       SHP.geometry = SHP.geometry.to_crs(CRS_Axes_proj4)

       Fig, Axes = plt.subplots(nrows=1, ncols=1, subplot_kw={'projection': CRS_Axes})

       SHP.plot(ax=Axes, transform=CRS_Axes)

End of Code

I thank you for your time, and I hope hearing from you soon.

Sincerely yours,

Philipe Riskalla Leal INSTITUTO NACIONAL DE PESQUISAS ESPACIAIS (INPE) - BRAZIL.

UofTCoders / studyGroup