Plotting geographic data using coordinates

Question

While it's possible to make a map plot by either:

crs = area_def.to_cartopy_crs()
ax = plt.axes(projection=crs)
ax.background_img(name='BM')  
plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper', cmap='RdBu_r')  

or

sst = dataset.variables['sst'][0, :, :]
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]

ax = plt.axes(projection=ccrs.PlateCarree())

plt.contourf(lons, lats, sst, 60,
             transform=ccrs.PlateCarree())

The first needs a crs object. And the second only works for filled contour plots. Is there a way get an imshow map plot with three arrays, data, lats, and lons.

Answer 1

Although you have cartopy in your tags, I think what you are trying to achieve can be solved with geopandas. The important concept is to have the same CRS as you plot your points in the figure so that all information align.

Lets look at a simple example

import geopandas
from matplotlib import pyplot as plt

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
cities = geopandas.read_file(geopandas.datasets.get_path('naturalearth_cities'))

ax = world.plot(color='white', edgecolor='black')
cities.plot(ax=ax, marker='o', color='red', markersize=5)
plt.show()

Note: Since we want to plot the cities on the same map we use the same figure Axes ax. Also note that both world and cities have the same CRS. You can see this by doing

print(world.crs, cities.crs)
epsg:4326 epsg:4326

Both return epsg:4326, so same CRS.

Now, you have a new set of points you want to add to your plot. Let's create a few random points.

from shapely import Point
import numpy as np

np.random.seed(1)
my_points = geopandas.GeoDataFrame(
    geometry=[Point(x, y) for x, y in zip(
        np.random.uniform(low=30, high=40, size=10),
        np.random.uniform(low=-30, high=-10, size=10)
    )], crs=world.crs
)

Here we create random points between lon [30, 40] east and lat [10, 30] south. Note that I am copying the crs of world since it is epsg:4326.

If it was something else, we would initialise my_points with crs='epsg:4326' and then translate my_points to world.crs as follows

my_points.to_crs(crs=world.crs, inplace=True)

Finally we can plot on the same Axes

my_points.plot(ax=ax, marker='s', color='g', markersize=10)

For more, visit this page

Answer 2

from:

https://coderzcolumn.com/tutorials/data-science/cartopy-basic-maps-scatter-map-bubble-map-and-connection-map

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

ax = fig.add_subplot(1,1,1, projection=crs.Robinson())

ax.set_global()

ax.add_feature(cfeature.COASTLINE, edgecolor="tomato")
ax.add_feature(cfeature.BORDERS, edgecolor="tomato")
ax.gridlines()

plt.scatter(x=starbucks_locations.Longitude, y=starbucks_locations.Latitude,
            color="dodgerblue",
            s=1,
            alpha=0.5,
            transform=crs.PlateCarree()) ## Important

plt.show()

I also use marker="s" and s=5 to have big square markers, this way its displayed as an image.

Answer 3

Not really, and it wouldn't make sense I think. When using imshow, you can plot any 2D array that you want, regardless of what projection it is, or whether it is actually geographic data. Combining it with other data would require specifying the correct extent as well, using the coordinates of the array you're plotting.

Since imshow requires a regular array (grid), it's not necessary to also specify the coordinate for every element, just the outer extent is sufficient. If this is all you want to do, you don't need to use Cartopy and just Matplotlib is good enough.

But as soon as you want to mix in other geographic elements, country borders, lat/lat gridlines, or other data in a different projection, you'll probably want to use Cartopy because it's very convenient for that.

The same is true when plotting an irregular grid, at which point methods like contourf or pcolormesh become useful. Those do require explicit coordinates for every element. And similar as with imshow those methods don't require Cartopy and a projection to be specified. That's just convenient when you want to mix in other geographic data, potentially in another projection.

You can always avoid the dependency on Cartopy (or Basemap etc) by making sure all data that you plot shares the same projection. But in most instances that's just not worth the effort, at least for me personally.