Search code examples
pythonmatplotlibcartopy

Plot square Cartopy map

I need to plot a square map using Cartopy. I currently use the following code for my map:

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

img = cimgt.GoogleTiles()

ax = plt.axes(projection = img.crs)
ax.set_extent((d['longitude'].min() - 0.05, d['longitude'].max() + 0.05,
               d['latitude'].min() - 0.05, d['latitude'].max() + 0.05))

ax.add_image(img, 10, interpolation = 'bicubic')

plt.scatter(d['longitude'], d['latitude'], transform = ccrs.PlateCarree(),
            c = '#E8175D', s = 14)

This works fine, except for the fact that the map isn't square. Instead, it's just fitted into the (15, 15) square of the plot.

enter image description here

I would like to add a bit more map to the left and to the right to make the plot perfectly square without distorting it. Simply setting the extent to the same difference on latitude and longitude doesn't do the job, because latitude and longitude cover different distances in Google's (and most other) map projections. I also found this post, but from what I get, the intent here is to distort the map.

I hope someone has an idea how to do this. It seems that Cartopy is not very intuitive in this regard.

Solution

  • To get square extent you need to specify it in map projection coordinates. That involves some coordinate transformation. Here is the code snippet that you need.

    # crs of your choice
crg = cimgt.StamenTerrain().crs    # or cimgt.GoogleTiles().crs

# set map limits, in degrees
lonmin, lonmax = -22, -15
latmin, latmax = 63, 65

# do coordinate transformation
LL = crg.transform_point(lonmin, latmin, ccrs.Geodetic())
UR = crg.transform_point(lonmax, latmax, ccrs.Geodetic())
EW = UR[0] - LL[0]
SN = UR[1] - LL[1]

# get side of the square extent (in map units, usually meters)
side = max(EW, SN)    # larger value is in effect
mid_x, mid_y = LL[0]+EW/2.0, LL[1]+SN/2.0  # center location

# the extent preserves the center location
extent = [mid_x-side/2.0, mid_x+side/2.0, mid_y-side/2.0, mid_y+side/2.0]

# this sets square extent
# crs=crg signifies that projection coordinates is used in extent
ax.set_extent(extent, crs=crg)

    Hope it helps.

    Edit

    Here is a complete working code and its resulting map.

    import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.io.img_tiles as cimgt

def make_map(projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(10, 10),
                           subplot_kw=dict(projection=projection))
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax

request = cimgt.StamenTerrain()   # very responsive
crg = request.crs   #crs of the projection
fig, ax = make_map(projection = crg)

# specify map extent here
lonmin, lonmax = -22, -15
latmin, latmax = 63, 65

LL = crg.transform_point(lonmin, latmin, ccrs.Geodetic())
UR = crg.transform_point(lonmax, latmax, ccrs.Geodetic())
EW = UR[0] - LL[0]
SN = UR[1] - LL[1]
side = max(EW,SN)
mid_x, mid_y = LL[0]+EW/2.0, LL[1]+SN/2.0  #center location

extent = [mid_x-side/2.0, mid_x+side/2.0, mid_y-side/2.0, mid_y+side/2.0]   # map coordinates, meters

ax.set_extent(extent, crs=crg)
ax.add_image(request, 8)

# add a marker at center of the map
plt.plot(mid_x, mid_y, marker='o', \
         color='red', markersize=10, \
         alpha=0.7, transform = crg)

plt.show()

    south of iceland