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