Getting contextily basemap to fill plots

Question

I'm trying to create a 2x2 subplot using the following:

import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx

site = gdf.groupby('Site_name')

plt.figure(figsize =(20,20))

# Iterate through sites

for i, (Site_name, site_gdf) in enumerate(site):
    # create subplot axes in a 2x2 grid
    ax = plt.subplot(2, 2, i + 1) # nrows, ncols, axes position
    # plot the site on these axes
    site_gdf.plot(ax=ax)
    ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery)
    # set the title
    ax.set_title(Site_name)
    # set the aspect
    # adjustable datalim ensure that the plots have the same axes size
    ax.set_aspect('equal', adjustable='datalim')

plt.tight_layout()
plt.show()

The problem occurs that the extent of the basemaps within the subplots are restricted to the overlaid point data. How do I change it so that each subplot has a basemap that covers the whole subplot.

Answer 1

I saw this post and a couple other similar ones in search for the same answer and came up with this as a way to take the user out of the loop in J.P.'s accepted answer.

This bit of code takes the desired figure size height and width, compares the orientation of your figure to the orientation of your gdf, and scales one of the axes "outward" to match the extent just as described in JP's answer above, just without having to adjust the scaling number manually. There is definitely some room for improvement (like defining a function or two), but I am looking forward to hearing about any holes you all find in this. It'll be a daily driver for me no doubt.

Best,

-Matt

import geopandas as gpd
import contextily as cx
import matplotlib.pyplot as plt

#Import Primary Geo Data
df = gpd.read_file(r'Primary.shp')

#import Secondary Layer
adf = gpd.read_file('Secondary.shp')
# - Reproject secondary layer if needed to match primary layer
if adf.crs != df.crs:
    adf = adf.to_crs(df.crs)
    print("Reprojected ADF")

#Plot Geodata
fig, ax = plt.subplots ()
width = 10
height = 15
if width > height:
    framedim = "wide"
else:
    framedim = "tall"
print(framedim)
plt.gcf().set_size_inches(width, height)
df.plot(ax = ax, alpha=0.5, edgecolor='silver', column = 'hour', cmap = "plasma", scheme=None, legend=False, zorder=3)

# - Set mapping extent to primary layer and user-defined figsize
lon_min, lat_min, lon_max, lat_max = df.total_bounds
dataWt = lon_max-lon_min
dataHt = lat_max-lat_min
print("dataWt: "+str(dataWt))
print("dataHt: "+str(dataHt))
if dataWt > dataHt:
    datadim = "wide"
else:
    datadim = "tall"
print(datadim)
#scaling
if (framedim == "tall") and (datadim == "tall"):
    #scale height
    print("Match Height")
    # - Get datawt unit for 1 unit of frame width
    one_data_unit = dataWt/width
    print("one_data_unit: "+str(one_data_unit))
    properheight = one_data_unit*height
    print("properheight: "+str(properheight))
    height_bookends = abs(properheight-dataHt)/2
    print("height_bookends: "+str(height_bookends))
    print("lon_min: "+str(lon_min))
    print("lon_max: "+str(lon_max))
    print("lat_min: "+str(lat_min))
    print("lat_max: "+str(lat_max))
    print("lat_min-height_bookends: "+str(lat_min-height_bookends))
    print("lat_max+height_bookends: "+str(lat_max+height_bookends))
    ax.set_xlim(lon_min-(one_data_unit*(width/100)), lon_max+(one_data_unit*(width/100)))
    ax.set_ylim(lat_min-height_bookends-(one_data_unit*(height/100)), lat_max+height_bookends+(one_data_unit*(height/100)))
elif (framedim == "wide") and (datadim == "wide"):
    #scale width
    print("Match Width")
    # - Get dataht unit for 1 unit of frame height
    one_data_unit = dataHt/height
    print("one_data_unit: "+str(one_data_unit))
    properwidth = one_data_unit*width
    print("properwidth: "+str(properwidth))
    width_bookends = abs(properwidth-dataWt)/2
    print("width_bookends: "+str(width_bookends))
    print("lon_min: "+str(lon_min))
    print("lon_max: "+str(lon_max))
    print("lon_min-width_bookends: "+str(lon_min-width_bookends))
    print("lon_max+width_bookends: "+str(lon_max+width_bookends))
    print("lat_min: "+str(lat_min))
    print("lat_max: "+str(lat_max))
    ax.set_xlim((lon_min-width_bookends)-(one_data_unit*(width/100)), (lon_max+width_bookends)+(one_data_unit*(width/100)))
    ax.set_ylim(lat_min-(one_data_unit*(height/100)), lat_max+(one_data_unit*(height/100)))
elif (framedim == "tall") and (datadim == "wide"):
    #scale height
    print("Fit Wide in Tall")
    # - Get datawt unit for 1 unit of frame width
    one_data_unit = dataWt/width
    print("one_data_unit: "+str(one_data_unit))
    properheight = one_data_unit*height
    print("properheight: "+str(properheight))
    height_bookends = abs(properheight-dataHt)/2
    print("height_bookends: "+str(height_bookends))
    print("lon_min: "+str(lon_min))
    print("lon_max: "+str(lon_max))
    print("lat_min: "+str(lat_min))
    print("lat_max: "+str(lat_max))
    print("lat_min-height_bookends: "+str(lat_min-height_bookends))
    print("lat_max+height_bookends: "+str(lat_max+height_bookends))
    ax.set_xlim(lon_min-(one_data_unit*(width/100)), lon_max+(one_data_unit*(width/100)))
    ax.set_ylim(lat_min-height_bookends-(one_data_unit*(height/100)), lat_max+height_bookends+(one_data_unit*(height/100)))
else:
    #scale width
    print("Fit Tall in Wide")
    # - Get dataht unit for 1 unit of frame height
    one_data_unit = dataHt/height
    print("one_data_unit: "+str(one_data_unit))
    properwidth = one_data_unit*width
    print("properwidth: "+str(properwidth))
    width_bookends = abs(properwidth-dataWt)/2
    print("width_bookends: "+str(width_bookends))
    print("lon_min: "+str(lon_min))
    print("lon_max: "+str(lon_max))
    print("lon_min-width_bookends: "+str(lon_min-width_bookends))
    print("lon_max+width_bookends: "+str(lon_max+width_bookends))
    print("lat_min: "+str(lat_min))
    print("lat_max: "+str(lat_max))
    ax.set_xlim((lon_min-width_bookends)-(one_data_unit*(width/100)), (lon_max+width_bookends)+(one_data_unit*(width/100)))
    ax.set_ylim(lat_min-(one_data_unit*(height/100)), lat_max+(one_data_unit*(height/100)))

# - Plot Secondary Layer now so that the pri layer is the focus of the extent
adf.plot(ax = ax, alpha=0.5, edgecolor='k', color = "r", zorder=2)

#Add in Basemap based on user definition
cx.add_basemap(ax, 
               crs=df.crs, 
               source=cx.providers.CartoDB.DarkMatterOnlyLabels, 
               attribution = False, 
               reset_extent=True
              )
cx.add_basemap(ax, 
               crs=df.crs, 
               source=cx.providers.CartoDB.DarkMatterNoLabels, 
               attribution = False,
               alpha = .8, 
               reset_extent=True
              )
ax.set_axis_off()

#Add in Marginalia
plt.text(.01, .03, str(df.crs).upper(), fontsize=10, color = "k", backgroundcolor = "gray", ha='left', va='top', transform=ax.transAxes, zorder=4)

plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
plt.margins(0,0)

fig.savefig("contextily.png")

Answer 2

my solution was, that i searched for my outermost points and then changed my matplotlib axes with

ax.set_xlim(x_min-addition, x_max+addition)
ax.set_ylim(y_min-addition, y_max+addition)

addition in this case is just a number to increase the view on the points, like 0.001. x_min was the smallest point, x_max the biggest, and so on... Afterwards, i got my basemap.

ctx.add_basemap(ax,...)

If you know your outermost boundaries for all points, you could pass them directly. The result should be in the same size. If you don't know your boundaries, you need to find them first. Therefore you could iterate over your points and just compare them (just an example):

for point_data_x, point_data_y in site_data:
    if point_data_x < x_min:
        x_min = point_data_x
    elif point_data_x > x_max:
        x_max = point_data_x
    if point_data_y < y_min:
        y_min = point_data_y
    elif point_data_y > y_max:
        y_max = point_data_y

Kind regards