feat: add Arctic/Greenland to utils/maps

src/polartoolkit/maps.py

@@ -1227,6 +1227,7 @@ def add_box(
 
 
 def interactive_map(
+    hemisphere: str,
     center_yx: list[float] | None = None,
     zoom: float = 0,
     display_xy: bool = True,
@@ -1237,10 +1238,12 @@ def interactive_map(
 ) -> ipyleaflet.Map:
     """
     Plot an interactive map with satellite imagery. Clicking gives the cursor location
-    in EPSG:3031 [x,y]. Requires ipyleaflet
+    in a Polar Stereographic projection [x,y]. Requires ipyleaflet
 
     Parameters
     ----------
+    hemisphere : str
+        choose between plotting in the "north" or "south" hemispheres
     center_yx : list, optional
         choose center coordinates in EPSG3031 [y,x], by default [0,0]
     zoom : float, optional
@@ -1283,7 +1286,11 @@ def interactive_map(
             center_ll = [points.lon.mean(), points.lat.mean()]
         else:
             # convert points to lat lon
-            points_ll: pd.DataFrame = utils.epsg3031_to_latlon(points)
+            if hemisphere == "south":
+                points_ll: pd.DataFrame = utils.epsg3031_to_latlon(points)
+            elif hemisphere == "north":
+                points_ll = utils.epsg3413_to_latlon(points)
+
             # if points supplied, center map on points
             center_ll = [np.nanmedian(points_ll.lat), np.nanmedian(points_ll.lon)]
             # add points to geodataframe
@@ -1293,25 +1300,47 @@ def interactive_map(
             )
             geo_data = ipyleaflet.GeoData(
                 geo_dataframe=gdf,
-                # style={'radius': .5, 'color': 'red', 'weight': .5},
                 point_style={"radius": 1, "color": "red", "weight": 1},
             )
     else:
         # if no points, center map on 0, 0
-        center_ll = utils.epsg3031_to_latlon([0, 0])
+        if hemisphere == "south":
+            center_ll = utils.epsg3031_to_latlon([0, 0])
+        elif hemisphere == "north":
+            center_ll = utils.epsg3413_to_latlon([0, 0])
 
     if center_yx is not None:
-        center_ll = utils.epsg3031_to_latlon(center_yx)
-
-    if basemap_type == "BlueMarble":
-        base = ipyleaflet.basemaps.NASAGIBS.BlueMarble3031  # pylint: disable=no-member
-        proj = ipyleaflet.projections.EPSG3031.NASAGIBS
-    elif basemap_type == "Imagery":
-        base = ipyleaflet.basemaps.Esri.AntarcticImagery  # pylint: disable=no-member
-        proj = ipyleaflet.projections.EPSG3031.ESRIImagery
-    elif basemap_type == "Basemap":
-        base = ipyleaflet.basemaps.Esri.AntarcticBasemap  # pylint: disable=no-member
-        proj = ipyleaflet.projections.EPSG3031.ESRIBasemap
+        if hemisphere == "south":
+            center_ll = utils.epsg3031_to_latlon(center_yx)
+        elif hemisphere == "north":
+            center_ll = utils.epsg3413_to_latlon(center_yx)
+
+    if hemisphere == "south":
+        if basemap_type == "BlueMarble":
+            base = ipyleaflet.basemaps.NASAGIBS.BlueMarble3031  # pylint: disable=no-member
+            proj = ipyleaflet.projections.EPSG3031.NASAGIBS
+        elif basemap_type == "Imagery":
+            base = ipyleaflet.basemaps.Esri.AntarcticImagery  # pylint: disable=no-member
+            proj = ipyleaflet.projections.EPSG3031.ESRIImagery
+        elif basemap_type == "Basemap":
+            base = ipyleaflet.basemaps.Esri.AntarcticBasemap  # pylint: disable=no-member
+            proj = ipyleaflet.projections.EPSG3031.ESRIBasemap
+        else:
+            msg = "invalid string for basemap_type"
+            raise ValueError(msg)
+    elif hemisphere == "north":
+        if basemap_type == "BlueMarble":
+            base = ipyleaflet.basemaps.NASAGIBS.BlueMarble3413  # pylint: disable=no-member
+            proj = ipyleaflet.projections.EPSG3413.NASAGIBS
+        # elif basemap_type == "Imagery":
+        #     base = ipyleaflet.basemaps.Esri.ArcticImagery  # pylint: disable=no-member
+        #     proj = ipyleaflet.projections.EPSG5936.ESRIImagery
+        # elif basemap_type == "Basemap":
+        #     base = ipyleaflet.basemaps.Esri.OceanBasemap  # pylint: disable=no-member
+        #     proj = ipyleaflet.projections.EPSG5936.ESRIBasemap
+        else:
+            msg = "invalid string for basemap_type"
+            raise ValueError(msg)
 
     # create the map
     m = ipyleaflet.Map(
@@ -1334,7 +1363,10 @@ def interactive_map(
         def handle_click(**kwargs: typing.Any) -> None:
             if kwargs.get("type") == "click":
                 latlon = kwargs.get("coordinates")
-                label_xy.value = str(utils.latlon_to_epsg3031(latlon))
+                if hemisphere == "south":
+                    label_xy.value = str(utils.latlon_to_epsg3031(latlon))
+                elif hemisphere == "north":
+                    label_xy.value = str(utils.latlon_to_epsg3413(latlon))
 
     m.on_interaction(handle_click)
 
@@ -1561,6 +1593,7 @@ def plot_3d(
             grid = utils.mask_from_polygon(
                 polygon_mask,
                 grid=grid,
+                hemisphere=kwargs.get("hemisphere", None),
             )
         # create colorscales
         if grd2cpt is True:

src/polartoolkit/utils.py

@@ -237,13 +237,17 @@ def region_to_df(
 
 
 def region_xy_to_ll(
-    region: tuple[typing.Any, typing.Any, typing.Any, typing.Any], dms: bool = False
+    region: tuple[typing.Any, typing.Any, typing.Any, typing.Any],
+    hemisphere: str,
+    dms: bool = False,
 ) -> tuple[typing.Any, typing.Any, typing.Any, typing.Any]:
     """
     Convert GMT region in format [e, w, n, s] in EPSG:3031 to lat / lon
 
     Parameters
     ----------
+    hemisphere : str,
+        choose between the "north" or "south" hemispheres
     region : tuple[typing.Any, typing.Any, typing.Any, typing.Any]
         region boundaries in GMT format; [e, w, n, s] in meters
     dms: bool, False
@@ -255,7 +259,10 @@ def region_xy_to_ll(
         region boundaries in GMT format; [e, w, n, s] in lat, lon
     """
     df = region_to_df(region)
-    df_proj = epsg3031_to_latlon(df, reg=True)
+    if hemisphere == "north":
+        df_proj = epsg3413_to_latlon(df, reg=True)
+    elif hemisphere == "south":
+        df_proj = epsg3031_to_latlon(df, reg=True)
 
     return tuple([dd2dms(x) for x in df_proj] if dms is True else df_proj)
 
@@ -391,6 +398,110 @@ def epsg3031_to_latlon(
     return df_out
 
 
+def latlon_to_epsg3413(
+    df: pd.DataFrame | NDArray[typing.Any, typing.Any],
+    reg: bool = False,
+    input_coord_names: tuple[str, str] = ("lon", "lat"),
+    output_coord_names: tuple[str, str] = ("x", "y"),
+) -> pd.DataFrame | NDArray[typing.Any, typing.Any]:
+    """
+    Convert coordinates from EPSG:4326 WGS84 in decimal degrees to EPSG:3413 North Polar
+    Stereographic in meters.
+
+    Parameters
+    ----------
+    df : pd.DataFrame or NDArray[typing.Any, typing.Any]
+        input dataframe with latitude and longitude columns
+    reg : bool, optional
+        if true, returns a GMT formatted region string, by default False
+    input_coord_names : list, optional
+        set names for input coordinate columns, by default ["lon", "lat"]
+    output_coord_names : list, optional
+        set names for output coordinate columns, by default ["x", "y"]
+
+    Returns
+    -------
+    pd.DataFrame or NDArray[typing.Any, typing.Any]
+        Updated dataframe with new easting and northing columns or NDArray in format
+        [e, w, n, s]
+    """
+    transformer = Transformer.from_crs("epsg:4326", "epsg:3413")
+
+    if isinstance(df, pd.DataFrame):
+        df_new = df.copy()
+        (  # pylint: disable=unpacking-non-sequence
+            df_new[output_coord_names[0]],
+            df_new[output_coord_names[1]],
+        ) = transformer.transform(
+            df_new[input_coord_names[1]].tolist(), df_new[input_coord_names[0]].tolist()
+        )
+    else:
+        ll = df.copy()
+        df_new = list(transformer.transform(ll[0], ll[1]))
+
+    if reg is True:
+        df_new = [
+            df_new[output_coord_names[0]].min(),
+            df_new[output_coord_names[0]].max(),
+            df_new[output_coord_names[1]].min(),
+            df_new[output_coord_names[1]].max(),
+        ]
+
+    return df_new
+
+
+def epsg3413_to_latlon(
+    df: pd.DataFrame | list[typing.Any],
+    reg: bool = False,
+    input_coord_names: tuple[str, str] = ("x", "y"),
+    output_coord_names: tuple[str, str] = ("lon", "lat"),
+) -> pd.DataFrame | list[typing.Any]:
+    """
+    Convert coordinates from EPSG:3413 North Polar Stereographic in meters to
+    EPSG:4326 WGS84 in decimal degrees.
+
+    Parameters
+    ----------
+    df : pd.DataFrame or list[typing.Any]
+        input dataframe with easting and northing columns, or list [x,y]
+    reg : bool, optional
+        if true, returns a GMT formatted region string, by default False
+    input_coord_names : list, optional
+        set names for input coordinate columns, by default ["x", "y"]
+    output_coord_names : list, optional
+        set names for output coordinate columns, by default ["lon", "lat"]
+
+    Returns
+    -------
+    pd.DataFrame or list[typing.Any]
+        Updated dataframe with new latitude and longitude columns, NDArray in
+        format [e, w, n, s], or list in format [lat, lon]
+    """
+
+    transformer = Transformer.from_crs("epsg:3413", "epsg:4326")
+
+    df_out = df.copy()
+
+    if isinstance(df, pd.DataFrame):
+        (  # pylint: disable=unpacking-non-sequence
+            df_out[output_coord_names[1]],  # type: ignore[call-overload]
+            df_out[output_coord_names[0]],  # type: ignore[call-overload]
+        ) = transformer.transform(
+            df_out[input_coord_names[0]].tolist(),  # type: ignore[call-overload]
+            df_out[input_coord_names[1]].tolist(),  # type: ignore[call-overload]
+        )
+        if reg is True:
+            df_out = [
+                df_out[output_coord_names[0]].min(),  # type: ignore[call-overload]
+                df_out[output_coord_names[0]].max(),  # type: ignore[call-overload]
+                df_out[output_coord_names[1]].min(),  # type: ignore[call-overload]
+                df_out[output_coord_names[1]].max(),  # type: ignore[call-overload]
+            ]
+    else:
+        df_out = list(transformer.transform(df_out[0], df_out[1]))
+    return df_out
+
+
 def points_inside_region(
     df: pd.DataFrame,
     region: tuple[float, float, float, float],
@@ -1543,13 +1654,18 @@ def get_min_max(
     return (v_min, v_max)
 
 
-def shapes_to_df(shapes: list[float]) -> pd.DataFrame:
+def shapes_to_df(
+    shapes: list[float],
+    hemisphere: str,
+) -> pd.DataFrame:
     """
     convert the output of `regions.draw_region` and `profile.draw_lines` to a dataframe
     of x and y points
 
     Parameters
     ----------
+    hemisphere : str
+        choose between the "north" or "south" hemispheres
     shapes : list
         list of vertices
 
@@ -1566,25 +1682,39 @@ def shapes_to_df(shapes: list[float]) -> pd.DataFrame:
         shape = pd.DataFrame({"lon": lon, "lat": lat, "shape_num": i})
         df = pd.concat((df, shape))
 
-    return latlon_to_epsg3031(df)
+    if hemisphere == "north":
+        df = latlon_to_epsg3413(df)
+    elif hemisphere == "south":
+        df = latlon_to_epsg3031(df)
+    else:
+        msg = "hemisphere must be 'north' or 'south'"
+        raise ValueError(msg)
+
+    return df
 
 
-def polygon_to_region(polygon: list[float]) -> tuple[float, float, float, float]:
+def polygon_to_region(
+    polygon: list[float],
+    hemisphere: str,
+) -> tuple[float, float, float, float]:
     """
-    convert the output of `regions.draw_region` to bounding region in EPSG:3031
+    convert the output of `regions.draw_region` to bounding region in EPSG:3031 for the
+    south hemisphere and EPSG:3413 for the north hemisphere.
 
     Parameters
     ----------
     polyon : list
         list of polygon vertices
+    hemisphere : str
+        choose between the "north" or "south" hemispheres
 
     Returns
     -------
     tuple[float, float, float, float]
         region in format [e,w,n,s]
     """
 
-    df = shapes_to_df(polygon)
+    df = shapes_to_df(shapes=polygon, hemisphere=hemisphere)
 
     if df.shape_num.max() > 0:
         logging.info(
@@ -1599,6 +1729,7 @@ def polygon_to_region(polygon: list[float]) -> tuple[float, float, float, float]
 
 def mask_from_polygon(
     polygon: list[float],
+    hemisphere: str,
     invert: bool = False,
     drop_nans: bool = False,
     grid: str | xr.DataArray | None = None,
@@ -1613,6 +1744,8 @@ def mask_from_polygon(
     ----------
     polygon : list
        list of polygon vertices
+    hemisphere : str
+        choose between the "north" or "south" hemispheres
     invert : bool, optional
         reverse the sense of masking, by default False
     drop_nans : bool, optional
@@ -1631,7 +1764,7 @@ def mask_from_polygon(
     """
 
     # convert drawn polygon into dataframe
-    df = shapes_to_df(polygon)
+    df = shapes_to_df(polygon, hemisphere=hemisphere)
     data_coords = (df.x, df.y)
 
     # remove additional polygons