-1

I'm migrating from basemap to cartopy. One thing I would like to do is plot high/low pressure on a map, such as in basemap. There is a good example on this page of how to do this: https://matplotlib.org/basemap/users/examples.html ("Plot sea-level pressure weather map with labelled highs and lows"). I'm not going to copy and paste the code from this site, but would like to know how to do the same in cartopy. The main thing I can't get my head around is how to do m.xmax and x > m.xmin and y < m.ymax and y > m.ymin in cartopy (some kind of vector transform I'd imagine.

I've had a good look and can't see this particular example translated into something compatible with cartopy. Any help would be welcome!

J W
  • 617
  • 1
  • 9
  • 28
  • The page you link to contains 10 examples. Are you asking for translating all of them to cartopy? – ImportanceOfBeingErnest Feb 27 '18 at 10:52
  • No, sorry, just the "Plot sea-level pressure weather map with labelled highs and lows" script – J W Feb 27 '18 at 10:55
  • The line in question is filtering out points that aren't in the projected map domain, so you could use the `get_extent()` method of a Cartopy axes to determine this domain (or you could make an initial implementation that doesn't bother with the filtering). – ajdawson Feb 27 '18 at 13:35
  • Thanks Andrew! Do you know if this example will appear in the Cartopy - as I know a few others who will want to know how to plot lows and highs on maps in the future using Cartopy.. http://scitools.org.uk/cartopy/docs/v0.15/gallery.html – J W Feb 27 '18 at 13:47
  • I'm not sure what the equivalent of x,y = m(x,y) is in Cartopy when converting from coordinates to grid space - it throws an error when I try transform_points() .. `x, y = np.meshgrid(LON, LAT)` ; `x, y = CRS.transform_points(ccrs.PlateCarree(),x, y)`, asking for 3 positional arguments (?) – J W Feb 27 '18 at 14:13
  • You are using the class, you need to call the `transform_points` method on an instance of a coordinate reference system. See my answer below. – ajdawson Feb 27 '18 at 14:37

1 Answers1

2

In order to write an equivalent program using cartopy you need to be able to translate two concepts. The first is finding the extent of a projection, this can be done with the get_extent() method of a GeoAxes:

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

my_proj = ccrs.Miller(central_longitude=180)
ax = plt.axes(projection=my_proj)
xmin, xmax, ymin, ymax = ax.get_extent()

You also need to transform coordinate points from geographic to projection coordinates, which is the function of the transform_points() method of a coordinate reference system instance:

import numpy as np

lons2d, lats2d = np.meshgrid(lons, lats) # lons lats are in degrees
transformed = my_proj.transform_points(ccrs.Geodetic(), lons2d, lats2d)
x = transformed[..., 0]  # lons in projection coordinates
y = transformed[..., 1]  # lats in projection coordinates

Now you can use the same technique as in the basemap example to filter and plot points, where instead of m.xmin you use xmin etc.

There are of course alternate ways of doing this which have pros and cons relative to the basemap example. If you come up with something nice you can contribute it to the Cartopy gallery.

ajdawson
  • 3,183
  • 27
  • 36
  • Great, thank you! Line `transformed = my_proj.transform_points(ccrs.Geodetic, lons2d, lats2d)` was what I was stuck on, now working. – J W Feb 27 '18 at 15:04