[FIXED] Struggling to plot data using cartopy and matplotlib

Issue

I am struggling with plotting points on top of a basemap using cartopy (and matplotlib). I think the problem is that I can't manage to get my data and the basemap to match projection/crs for the axes, as my data structure is different from most examples I found. Unforunatly some of it (top part in my code) can't be changed, as it would mean rewriting a lot of other code and functions :/

At first I tried using matplotlib to plot the points and cartopy to add the basemap (Version 1 in my code). Depending on where I set the projection for ax the basemap is not drawn at all or on top of everything without the wanted extend (-> whole earth is shown). If I add the transform= parameter as suggested here, nothing gets drawn at all (at least no inside the viewport).

On geopandas.org I found an example where a cartopy-function (add_geometries) is used to add the points (Version 2). This throws the Error Source CRS must be an instance of CRS or one of its subclasses, or None. where I don't know where it comes from.

Here is the code I use. Unfortunatly it's rather lengthy, but I didn't know what else I could omit without loosing necessary information.

import pandas as pd                            #2.1.4
import geopandas as gpd                        #0.12.2

import matplotlib.pyplot as plt                #3.8.0
from shapely.geometry import Point, Polygon    #1.8.5

import cartopy.crs as ccrs                     #0.22.0
import cartopy.io.img_tiles as cimgt

#------------------------------------------
#-----------Creating data-environment------
#------------------------------------------

#------Part 1: has to stay--------
gdf_test1 = gpd.GeoDataFrame(
    {
        "geometry": [Point(13.10157, 52.36640), Point(13.10130, 52.36656), Point(13.10254, 52.36582), Point(13.10309, 52.36549), Point(13.10332, 52.36535)]
    },
    crs = 'EPSG:4326'
)

gdf_test2 = gpd.GeoDataFrame(
    {
        "geometry": [Point(13.43562, 52.48038), Point(13.50074, 52.45651), Point(13.50061, 52.45657), Point(13.49998, 52.45688), Point(13.49984, 52.45695)]
    },
    crs = 'EPSG:4326'
)

outline_test1 = Polygon((Point(12.984381, 52.354036), Point(13.147663, 52.354036), Point(13.147663, 52.412748), Point(12.984381, 52.412748), Point(12.984381, 52.354036)))
outline_test2 = Polygon((Point(13.62631009242705, 52.517735730410585) , Point(13.625180130891495, 52.50452520302402) , Point(13.621934955922162, 52.49144577310617) , Point(13.616607639378886, 52.478623231164505) , Point(13.609251157772427, 52.46618078976451) , Point(13.599937818411506, 52.45423790766038) , Point(13.5887585076519, 52.44290915361389) , Point(13.57582177055312, 52.432303120286605) , Point(13.561252732467345, 52.42252139799845) , Point(13.545191874144964, 52.413657617485086) , Point(13.527793672849741, 52.40579657006711) , Point(13.509225122743361, 52.399013412873614) , Point(13.489664148436754, 52.39337296594993) , Point(13.469297926125792, 52.388929107231334) , Point(13.448321127145887, 52.38572427048734) , Point(13.426934099105035, 52.3837890504421) , Point(13.405341000000002, 52.38314191835826) , Point(13.383747900894965, 52.3837890504421) , Point(13.362360872854113, 52.38572427048734) , Point(13.341384073874208, 52.388929107231334) , Point(13.321017851563246, 52.39337296594993) , Point(13.301456877256639, 52.399013412873614) , Point(13.282888327150259, 52.40579657006711) , Point(13.265490125855036, 52.413657617485086) , Point(13.249429267532655, 52.42252139799845) , Point(13.23486022944688, 52.432303120286605) , Point(13.221923492348102, 52.44290915361389) , Point(13.210744181588494, 52.45423790766038) , Point(13.201430842227573, 52.46618078976451) , Point(13.194074360621114, 52.478623231164505) , Point(13.188747044077843, 52.49144577310617) , Point(13.185501869108505, 52.50452520302403) , Point(13.18437190757295, 52.517735730410585) , Point(13.185369938647879, 52.53095019147478) , Point(13.188488252991643, 52.5440412712636) , Point(13.193698653775746, 52.5568827315916) , Point(13.200952657395336, 52.569350632900345) , Point(13.210181894677728, 52.581324538061054) , Point(13.221298711298406, 52.59268868614871) , Point(13.234196963913654, 52.60333312435861) , Point(13.248753006259824, 52.61315478651281) , Point(13.264826857184044, 52.62205850701603) , Point(13.282263540302736, 52.62995795966922) , Point(13.300894582770653, 52.63677651143157) , Point(13.320539658531636, 52.6424479820343) , Point(13.341008359453616, 52.64691730128411) , Point(13.362102075970094, 52.65014105694213) , Point(13.383615967296576, 52.65208792721094) , Point(13.405341000000002, 52.652738993096214) , Point(13.427066032703426, 52.65208792721094) , Point(13.448579924029907, 52.65014105694215) , Point(13.469673640546384, 52.64691730128411) , Point(13.490142341468363, 52.642447982034305) , Point(13.509787417229353, 52.63677651143158) , Point(13.528418459697264, 52.62995795966922) , Point(13.54585514281595, 52.62205850701603) , Point(13.56192899374018, 52.61315478651281) , Point(13.57648503608635, 52.60333312435861) , Point(13.589383288701594, 52.592688686148705) , Point(13.600500105322267, 52.581324538061054) , Point(13.609729342604664, 52.56935063290036) , Point(13.616983346224249, 52.5568827315916) , Point(13.622193747008357, 52.5440412712636) , Point(13.625312061352123, 52.53095019147478) , Point(13.62631009242705, 52.517735730410585)))

liste_test = pd.DataFrame({
    'order': [1, 2],
    'bez': ['Region 1 (circ)', 'Region 2 (rect)'],
    'gdf': [gdf_test1, gdf_test2],
    'filter': [outline_test1, outline_test2]
})

#------Part 2: could be changed--------

gdf_lang_test = gpd.GeoDataFrame(columns=['bez', 'geometry'], geometry='geometry', crs='EPSG:4326')  
for i, row in liste_test.iterrows():
    for j, row2 in row['gdf'].iterrows():
        ptX = row2['geometry'].x
        ptY = row2['geometry'].y
        gdf_this_test = gpd.GeoDataFrame([(row['bez'], Point(ptX, ptY))], columns=['bez', 'geometry'], geometry='geometry', crs='EPSG:4326')        #jeder Messwert innerhalb eines gdf
        gdf_lang_test = gpd.GeoDataFrame(pd.concat([gdf_lang_test, gdf_this_test], ignore_index=True))  #wird an neuen gdf angehängt

bland_test1 = gpd.read_file('https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/2_bundeslaender/1_sehr_hoch.geo.json')
bland_test2 = gpd.read_file('https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/4_kreise/1_sehr_hoch.geo.json')

args = [bland_test1, bland_test2, gdf_lang_test]

#-----------------------------------------------
#-------------Actual Code-----------------------
#-----------------------------------------------

fig, ax = plt.subplots(figsize=(15, 15))
bland_crs = args[0].crs   #Returns: epsg:4326  (lower case!?)
data_crs = args[2].crs    #Returns: EPSG:4326  (upper case!?)

ax = plt.axes(projection=ccrs.Mercator())  

#-----PLOTTING DATA-----------------------------
#-Version 1----------
args[0].plot(ax=ax,facecolor="none", edgecolor='black', linewidths=1, alpha=1, figsize=(15,15))      #, transform=ccrs.PlateCarree())
args[1].plot(ax=ax,facecolor="none", edgecolor='grey', linewidths=0.5, alpha=1, figsize=(15,15))     #, transform=ccrs.PlateCarree()(
args[2].plot(ax=ax, zorder=2, markersize=6)     #, transform=ccrs.PlateCarree())

#-Version 2----------
#ax.add_geometries(args[0]['geometry'], facecolor="none", edgecolor='black', crs=bland_crs)
#ax.add_geometries(args[0]['geometry'], facecolor="none", edgecolor='black', crs=bland_crs)
#ax.add_geometries(args[2]['geometry'], facecolor="none", edgecolor='black', crs=data_crs)

#-----PLOTTING Outlines-------------------------
n=0
for i, row in liste_test.iterrows():
    n=n+1
    polygon_temp = row['filter']
    p = gpd.GeoSeries(polygon_temp, crs=data_crs)      
    args += (p,)                             
    #-Version 1----
    args[n+2].plot(ax=ax, facecolor="none")    #,transform=ccrs.PlateCarree())

    #-Version 2---
    #ax.add_geometries(args[n+2], facecolor="none", edgecolor='black', crs=data_crs)

#----Calculate Extent---------------------------
xmin, ymin, xmax, ymax = args[2].total_bounds
xmin = xmin*0.998
xmax = xmax*1.002
ymin = ymin*0.998
ymax = ymax*1.002
xlim = ([xmin, xmax])
ylim = ([ymin, ymax])
ax.set_xlim(xlim)
ax.set_ylim(ylim)

#----PLOTTING Background image--------------
#ax = plt.axes(projection=ccrs.Mercator())     #if uncommented, the basemap is drawn on top of everything with max extent (whole earth)

url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}.jpg'
ax.set_extent(ax.get_extent(), ccrs.Mercator())
image = cimgt.GoogleTiles(url=url)
ax.add_image(image, 1)

plt.show()

Help would be much appreciated! :/

Solution

You should probably take a step back and start simplifying your code. It already starts with the issue of creating two axes, one with and one without a Cartopy projection. Overwriting the ax variable doesn't make the fist axes disappear from the result.

I usually also avoid mixing two interfaces for plotting, while technically fine, it easily gets confusing about which library does what. In this case you use both Matplotlib and Geopandas, it's probably clearer to pick one. Either do all plotting with Matplotib (more explicit/low level) or all with Geopandas.

In this case using Matplotlib (like your 2nd version) has some performance benefits, I suspect because Cartopy is a bit smarter when reprojecting the detailed vectors to the tiny final extent. Most of that data doesn't end up on the plot, so for performance it might make sense to intersect it upfront.

And make sure to always assign the appropriate projection (transform= or crs=) when plotting. It might help to temporarily draw the gridlines/labels, and add elements one at a time, since that will quickly show where it fails.

The example below manually assigns the projection. Converting a Geopandas crs to Cartopy isn't that straightforward (I think...), but an automatic approach can be done by first converting the GeoDataframe to the map projection as shown at:
https://geopandas.org/en/stable/gallery/cartopy_convert.html#Plotting-with-CartoPy

fig, ax = plt.subplots(figsize=(6, 8), subplot_kw=dict(projection=ccrs.Mercator()))

xmin, ymin, xmax, ymax = gdf_lang_test.total_bounds
ax.set_extent(
    (xmin*0.998, xmax*1.002, ymin*0.998, ymax*1.002), 
    ccrs.PlateCarree(),
)

bg_image = cimgt.GoogleTiles(
    url='https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}.jpg',
)
ax.add_image(bg_image, 13, interpolation='hanning', regrid_shape=3000)

ax.add_geometries(bland_test1['geometry'], fc="none", ec='k', crs=ccrs.PlateCarree())
ax.add_geometries(bland_test2['geometry'], fc="none", ec='k', crs=ccrs.PlateCarree())
ax.plot(*zip(*[(p.x,p.y) for p in gdf_lang_test['geometry']]), "o", mec="k", mew=1, transform=ccrs.PlateCarree())

for i, row in liste_test.iterrows():    
    ax.add_geometries(row['filter'], fc="none", ec="r", crs=ccrs.PlateCarree())

ax.gridlines(draw_labels=True)

Answered By - Rutger Kassies