CODEWITHSUNDEEP
pythonmatplotlibcartopy
Will Mc.
Will Mc.

Reputation: 31

Issue w/ image crossing dateline in imshow & cartopy

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.

Upvotes: 3

Views: 1354

Answers (1)

swatchai
swatchai

Reputation: 18782

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

Upvotes: 1

Related Questions