Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Projection Problems when Displaying an Image on a Map with Cartopy

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 wrong projection

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.

no coastlines

I saw the issue reported here. But

ax.set_global()

results in this image:

where is my data

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:

correct image

like image 422
Julien Chastang Avatar asked Sep 10 '13 22:09

Julien Chastang


1 Answers

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

like image 120
pelson Avatar answered Jan 03 '23 13:01

pelson