Question

I received the following email and wanted to make sure the answer to this question was available to everybody:


Hi,

I would like to setup a simple latitude longitude map, using cartopy, which crosses the dateline and shows east Asia on the left hand side with the west of North America on the right. The following google map is roughly what I am after:

https://maps.google.co.uk/?ll=56.559482,-175.253906&spn=47.333523,133.066406&t=m&z=4

Screenshot of google map

Can this be done with Cartopy?

Was it helpful?

Solution

Good question. This is probably something which will come up time-and-time again, so I will go through this step-by-step before actually answering your specific question. For future reference, the following examples were written with cartopy v0.5.

Firstly, it is important to note that the default "latitude longitude" (or more technically PlateCarree) projection works in the forward range of -180 to 180. This means that you cannot plot the standard PlateCarree projection beyond this. There are several good reasons for this, most of which boil down to the fact that cartopy would have to do a lot more work when projecting both vectors and rasters (simple coastlines for example). Unfortunately the plot you are trying to produce requires precisely this functionality. To put this limitation into pictures, the default PlateCarree projection looks like:

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

proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.axes(projection=proj)
ax1.stock_img()
plt.title('Global')

plt.show()

global plate carree

Any single rectangle that you can draw on this map can legally be a zoomed in area (there is some slightly more advanced code in here, but the picture is worth a 1000 words):

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=-90, maxx=45, miny=15, maxy=70)
x0, y0, x1, y1 = box.bounds

proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.subplot(211, projection=proj)
ax1.stock_img()
ax1.add_geometries([box], proj, facecolor='coral', 
                   edgecolor='black', alpha=0.5)
plt.title('Global')

ax2 = plt.subplot(212, projection=proj)
ax2.stock_img()
ax2.set_extent([x0, x1, y0, y1], proj)
plt.title('Zoomed in area')

plt.show()

global with a second zoomed in map

Unfortunately the plot you want would require 2 rectangles with this projection:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=120, maxx=260, miny=15, maxy=80)

proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.axes(projection=proj)
ax1.stock_img()
ax1.add_geometries([box], proj, facecolor='coral', 
                   edgecolor='black', alpha=0.5)
plt.title('Target area')

plt.show()

target rectangle on global plot

Hence it is not possible to draw a map that crosses the dateline using the standard PlateCarree definition. Instead we could change the PlateCarree definition's central longitude to allow a single box to be drawn of the area we are targeting:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=120, maxx=260, miny=15, maxy=80)
x0, y0, x1, y1 = box.bounds

proj = ccrs.PlateCarree(central_longitude=180)
box_proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.subplot(211, projection=proj)
ax1.stock_img()
ax1.add_geometries([box], box_proj, facecolor='coral', 
                   edgecolor='black', alpha=0.5)
plt.title('Global')

ax2 = plt.subplot(212, projection=proj)
ax2.stock_img()
ax2.set_extent([x0, x1, y0, y1], box_proj)
plt.title('Zoomed in area')

plt.show()

Two PlateCarre plots with a central longitude of 180, one global, one zoomed in to the target area

Hopefully that shows you what it is you have to do to achieve your target map, the code above might be a little complex to achieve your goal, so to simplify slightly, the code I would write to produce the plot you want would be something like:

import cartopy.feature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt    

ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([120, 260, 15, 80], crs=ccrs.PlateCarree())

# add some features to make the map a little more polished
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines('50m')

plt.show()

final target map

This was a long answer, hopefully I have not only answered the question, but made some of the more complex details of map production and cartopy more clear to help smooth any future problems you may have.

Cheers,

OTHER TIPS

For more details of above benjimin's comment,

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.ticker import AutoMinorLocator, FixedLocator, MultipleLocator

def map_common(ax1,gl_loc=[True,True,False,True],gl_lon_info=range(-180,180,60),gl_dlat=30):

    ax1.coastlines(color='silver',linewidth=1.)

    gl = ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                       linewidth=0.6, color='gray', alpha=0.5, linestyle='--')

    gl.ylabels_left = gl_loc[0]
    gl.ylabels_right = gl_loc[1]
    gl.xlabels_top = gl_loc[2]
    gl.xlabels_bottom = gl_loc[3]

    gl.xlocator = FixedLocator(gl_lon_info)
    gl.ylocator = MultipleLocator(gl_dlat)
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 11, 'color': 'k'}
    gl.ylabel_style = {'size': 11, 'color': 'k'}


lon_boundary=np.arange(-240,-60,1.)
lat_boundary=np.arange(15,75,1.)
data=np.ones([lat_boundary.shape[0]-1,lon_boundary.shape[0]-1]) ## Data dimension is 1 less than boundaries
data=data*lat_boundary[:-1,None]

lon_offset=-150  ##
x,y=np.meshgrid(lon_boundary-lon_offset,lat_boundary)

fig=plt.figure()
fig.set_size_inches(7.5,5) ## (xsize, ysize)
ax1=fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=lon_offset))
ax1.set_extent([-250,-50,10,80],crs=ccrs.PlateCarree())

props=dict(vmin=0,vmax=90,cmap=plt.cm.get_cmap('bone'),alpha=0.8)
cs=ax1.pcolormesh(x,y,data,**props)
ax1.set_title('Lon_Offset=-90')
map_common(ax1,gl_lon_info=[-180,-120,-60,120,],gl_dlat=15)

fnout='./map_over_dateline.png'
#plt.show()
plt.savefig(fnout,bbox_inches='tight',dpi=150)

Output of this program

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top