Issue w/ image crossing dateline in imshow & cartopy

user_8952534 0
Viewed 1035

I'm trying to plot a square grid of equally-spaced (in lat/lon) data using cartopy, matplotlib, and imshow. The data crosses the dateline, and I've had issues getting a map to work properly.

Here's an example of my issue:

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


lat = np.arange(6000)*0.02 + (-59.99)
lon = np.arange(6000)*0.02 + (85.01)

dat = np.reshape(np.arange(6000*6000),[6000,6000])

tran = ccrs.PlateCarree()
proj = tran

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

ax = plt.axes(projection=proj)
print([lon[0],lon[-1],lat[0],lat[-1]])
ax.imshow(dat, extent=[lon[0],lon[-1],lat[0],lat[-1]],transform=tran,interpolation='nearest')
ax.coastlines(resolution='50m', color='black', linewidth=2)
ax.gridlines(crs=proj,draw_labels=True)
plt.show()

tran = ccrs.PlateCarree(central_longitude=180)
proj = tran

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

ax = plt.axes(projection=proj)
print([lon[0]-180,lon[-1]-180,lat[0],lat[-1]])
ax.imshow(dat, extent=[lon[0]-180,lon[-1]-180,lat[0],lat[-1]],transform=tran,interpolation='nearest')
ax.coastlines(resolution='50m', color='black', linewidth=2)
ax.gridlines(crs=tran,draw_labels=True)
plt.show()

The first plot yields this image, chopping off at 180E: chopped at 180

The second fixes the map issue, but the grid ticks are now wrong: 180 now zero

I've tried reprojecting, I think (where tran != proj), but it seemingly either hung or was taking too long.

I basically want the bottom image, but with the proper labels. I'm going to have more geolocated data to overplot, so I'd like to do it correctly instead of what seems like a hack right now.

cartopy matplotlib python
1 Answers

With Cartopy, drawing a map crossing dateline is always challenging. Here is the code that plots the map you want.

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

# demo raster data
n1 = 300
m1 = 0.4
lat = np.arange(n1)*m1 + (-59.99)
lon = np.arange(n1)*m1 + (85.01)
dat = np.reshape(np.arange(n1*n1), [n1,n1])

cm_lon=180  # for central meridian

tran = ccrs.PlateCarree(central_longitude = cm_lon)
proj = tran

plt.figure(figsize=(8,8))
ax = plt.axes(projection=proj)

ext = [lon[0]-cm_lon, lon[-1]-cm_lon, lat[0], lat[-1]]
#print(ext)

ax.imshow(dat, extent=ext, \
              transform=tran, interpolation='nearest')

ax.coastlines(resolution='110m', color='black', linewidth=0.5, zorder=10)

# this draws grid lines only, must go beyond E/W extents
ax.gridlines(draw_labels=False, xlocs=[80,100,120,140,160,180,-180,-160,-140])

# this draw lables only, exclude those outside E/W extents
ax.gridlines(draw_labels=True, xlocs=[100,120,140,160,180,-160])

plt.show()

The resulting map:

output

user_21774130
Related