10. Reproject your data

PolarToolkit works exclusively with data projected in either EPSG 3031 or EPSG 3413, for Antarctica and the Arctic/Greenland, respectively. These projections have units of meters and are relative to the poles. If you have your own data which you want to use in PolarToolkit, you may need to reproject it into these projections if its not already in them. Here we’ll demonstrate the basics of how to do this.

%load_ext autoreload
%autoreload 2
import os

from polartoolkit import fetch, maps, regions, utils

# set default to southern hemisphere for this notebook
os.environ["POLARTOOLKIT_HEMISPHERE"] = "south"

10.1. Convert a region from lat lon to EPSG

If you have a geographic region in latitude and longitude in the format (lon_min, lon_max, lat_min, lat_max) you can use region_ll_to_xy to project it to EPSG polar stereographic, or region_xy_to_ll for the opposite.

region_lat_lon = (-104.40, -92.33, -76.42, -73.89)
region_xy = utils.region_ll_to_xy(region_lat_lon)
region_xy

(np.float64(-1760092.1034694898),
 np.float64(-1435602.3038731467),
 np.float64(-438079.2907721378),
 np.float64(-60257.458812579534))

figure_region = regions.alter_region(region_xy, zoom=-200e3)

fig = maps.basemap(
    modis_basemap=True,
    region=figure_region,
)

fig.plot(
    x=[region_xy[0], region_xy[0], region_xy[1], region_xy[1], region_xy[0]],
    y=[region_xy[2], region_xy[3], region_xy[3], region_xy[2], region_xy[2]],
    pen="4p,white",
)


fig.plot(
    projection=utils.set_proj(figure_region, fig_height=15)[1],
    region=(*figure_region, "+ue", ""),  # codespell:ignore ue
    x=[
        region_lat_lon[0],
        region_lat_lon[0],
        region_lat_lon[1],
        region_lat_lon[1],
        region_lat_lon[0],
    ],
    y=[
        region_lat_lon[2],
        region_lat_lon[3],
        region_lat_lon[3],
        region_lat_lon[2],
        region_lat_lon[2],
    ],
    pen="2p,red",
)

fig.show(dpi=200)

Warning 1: The definition of geographic CRS EPSG:4326 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
Warning 1: The definition of geographic CRS EPSG:4326 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
Warning 1: The definition of geographic CRS EPSG:4326 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
Warning 1: The definition of geographic CRS EPSG:4326 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
Warning 1: The definition of geographic CRS EPSG:4326 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

The same can be done for coordinate pairs, or pandas dataframes of coordinates.

xy_coords = (-360000, -1100000)

lat_lon_coords = utils.epsg3031_to_latlon(xy_coords)

print(lat_lon_coords)

(-161.87813975209863, -79.37689147851952)

# plot on map
fig = maps.basemap(
    imagery_basemap=True,
    region=regions.roosevelt_island,
    inset=True,
)

fig.plot(
    x=xy_coords[0],
    y=xy_coords[1],
    style="c.6c",
    fill="black",
)
fig.plot(
    x=lat_lon_coords[0],
    y=lat_lon_coords[1],
    region=(
        *regions.roosevelt_island,
        "+ue",  # codespell:ignore ue
        "",
    ),  # inform region is in meters
    projection=utils.set_proj(
        regions.roosevelt_island,
        fig_height=15,
    )[1],
    style="t.8c",
    pen="2p,red",
)
fig.show(dpi=200)

If you have gridded data which needs to be reprojected, we recommend the rioxarray package with an example shown below.

# fetch a grid which is already in EPSG:3031
grid = fetch.sediment_thickness(
    version="GlobSed",
    spacing=5e3,
    region=regions.roosevelt_island,
    verbose="e",
)

# tell rioxarray what the current projection is
grid = grid.rio.write_crs("EPSG:3031")

# reproject to lat lon
grid_ll = grid.rio.reproject("EPSG:4326", resolution=0.05)

# reproject back to EPSG to see if it worked
grid_xy = grid_ll.rio.reproject("EPSG:3031", resolution=5e3)

_ = utils.grd_compare(
    grid,
    grid_xy,
    hist=True,
    robust=True,
    grid1_name="EPSG:3031",
    grid2_name="reprojected EPSG:3031",
)