I have some satellite image data I would like to display using Cartopy. I have successfully followed the image example detailed here. Resulting in this code:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig = plt.figure(figsize=(12, 12))
img_extent = (-77, -59, 9, 26)
ax = plt.axes(projection=ccrs.PlateCarree())
# image data coming from server, code not shown
ax.imshow(img, origin='upper', extent=img_extent)
ax.set_xmargin(0.05)
ax.set_ymargin(0.10)
# mark a known place to help us geo-locate ourselves
ax.plot(-117.1625, 32.715, 'bo', markersize=7)
ax.text(-117, 33, 'San Diego')
ax.coastlines()
ax.gridlines()
plt.show()
This code generates the following image
My problem is that the satellite image data is not in the PlateCarree projection, but the Mercator projection.
But when I get the axis object with
ax = plt.axes(projection=ccrs.Mercator())
I lose the coastlines.
I saw the issue reported here. But
ax.set_global()
results in this image:
The data is not present, and San Diego is in the wrong location. Also the lat/lon extents have changed. What am I doing wrong?
Post Discussion Update
The main problem is that I had not properly specified the image extents in the target projection with the transform_points
method. I also had to be specific about the coordinate reference system in the imshow
method as Phil suggests. Here is the correct code:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
proj = ccrs.Mercator()
fig = plt.figure(figsize=(12, 12))
extents = proj.transform_points(ccrs.Geodetic(),
np.array([-77, -59]),
np.array([9, 26]))
img_extents = (extents[0][0], extents[1][0], extents[0][6], extents[1][7] )
ax = plt.axes(projection=proj)
# image data coming from server, code not shown
ax.imshow(img, origin='upper', extent=img_extents,transform=proj)
ax.set_xmargin(0.05)
ax.set_ymargin(0.10)
# mark a known place to help us geo-locate ourselves
ax.plot(-117.1625, 32.715, 'bo', markersize=7, transform=ccrs.Geodetic())
ax.text(-117, 33, 'San Diego', transform=ccrs.Geodetic())
ax.coastlines()
ax.gridlines()
plt.show()
Resulting in this correctly geoprojected satellite image:
Ideally, try to always be specific about the coordinate reference system your data are in when plotting with cartopy (via the transform keyword). This will mean you can just switch projections in your script and the data will automatically be put in the correct place.
So in your case, the plt.imshow
should have a transform=ccrs.Mercator()
keyword argument (you may need a a more specific parameterised Mercator instance). If your extents are in Geodetic (lats and lons) you will have to transform the bounding box into the mercator coordinates, but other than that, everything else should work as expected.
NOTE: I'm going to go and update the example to include the transform argument ;-) (PR: https://github.com/SciTools/cartopy/pull/343)
HTH
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With