Geopandas: Earth as a Sphere

This one kept me busy for actually a few months. Map projections, the art of mapping the surface of a sphere onto a two dimensional canvas, is a science in its own rights. Dealing with data from the arctic will cause grief, as most maps we commonly use, limit the range of latitudes displayed to well below +/-90° latitude.

As an example, https://epsg.io/3857, WGS84, Spherical Mercator, has bounds of -85.06° and 85.06°. Not that many people live in the pole regions, after all.

Which is not the reason, the reason is the difficulty in projecting surfaces from the pole regions onto a two dimensional surface, the point at the poles would become a line. Which is why there is a range of projections defined for polar zones, (see Coordinate reference systems for “arctic” (epsg.io)).

For our mini apps on https://github.com/arctic-risk/arctic-data/tree/main/miniapps, we, however, need to combine geodata from various sources:

  • The amazing blue marble sea ice extent data set e.g. from NSIDC, and
  • the wonderfully curated earth map from naturalearth.com

Both come in different projections. If you, say, download the arctic sea ice extension shapefile from the http mirror of the NSIDC ftp site for September, 2020, then unzip and load it using the fabulous geopandas package, it will have EPSG:3411, NSIDC’s Sea Ice Polar Stereographic North, projection.

def load_ice_extent_shapefile(year,month):
    refdate = pd.to_datetime("{}-{}-15".format(wYear.value,wMonth.value))
    url = "http://masie_web.apps.nsidc.org/pub/DATASETS/NOAA/G02135/north/monthly/shapefiles/shp_extent/{:%m}_{:%b}/extent_N_{:%Y%m}_polygon_v3.0.zip".format(refdate,refdate,refdate)
    #return url
    DIR = "./sun"+urllib.parse.urlparse(url).path.replace(".zip","")
    if not os.path.exists(DIR):
        r = requests.get(url)
        if r.ok:
            zf = zipfile.ZipFile(io.BytesIO(r.content))
        Path(DIR).mkdir(exist_ok=True,parents=True)
        zf.extractall(DIR)
        
    for f in os.listdir(DIR):
        if ".shp" in f:
            gdfIceExtent = gpd.read_file(os.path.join(DIR,f))
    
    return gdfIceExtent

gdfArcticSeaIce202009 = load_ice_extent_shapefile(2020,9)
print(gdfArcticSeaIce202009.crs)
epsg:3411

Plotting it shows a shape that needs getting used to as it is a pancake style view of the north pole/arctic area, as this is sea ice extent, Greenland is missing, which gives a hint to what the map actually shows.

If we want to compare region sizes, we need a world map, one of the best sources of which being https://www.naturalearthdata.com/ . If you download and unzip the map data of your preference (I recommend 10m admin subunits from the cultural vectors download page), you can use these data to plot the world.

You can already tell from the axes that these two plots cannot be easily combined. You can plot one inside the other, but the arctic axes range are in the order of 1e6, while the world map range is the common +-180/+-90(ish).

So, we need to change the projection of the two maps, and align them. You can query the coordinate reference systems of the two maps by looking at the crs attribute of the GeoDataFrame object.

Now, you can convert GeoDataFrames’ projections by calling the .to_crs() method. If you do this, you will be surprised by the result:

At least that result makes slightly more sense than trying the other way around

The allowed range of coordinates do matter. So, we need a different projection. It took a few iterations until I came across CRS proj4 Strings, which allows for specifying an almost arbitrary projection. The string

+proj=ortho +lat_0=90 +lon_0=0 +x_0=0 +y_0=0

generates a “google Earth” style projection of a spherical earth seen from “the top”.

This finally resolved the issue around comparing the sea ice extent with regions and countries. The algorithm is quite straightforward:-

  • load the sea ice extent data
  • load the world data
  • apply the lat_0=90 lon_0=0 projection to the sea ice map
  • from the selected regions, compute the centroid of the map
  • use the map centroids as values for lat_0 and lon_0 of the region map and apply that projection. This way you “rotate the globe such that the centroid of the selected countries lines up with the north pole”

This is how the map overlay of the Ice Cover Comparisons mini-app is computed. The computation of the world (subregion) map centroid is not quite accurate and could be improved, but another improvement for that could be to add slide widgets to the streamlit app to allow for a user interaction with the map. After all, this mini app is a source of inspiration, so, why not leave room for improvement.