Fix projection in area plots (#1170)

* use mercator projection for standard area plots * include projection in area plot query * fix projections for polar plots
DFO-Ocean-Navigator · Nov 8, 2024 · f42b55a · f42b55a
1 parent aa5cacb
commit f42b55a
Show file tree

Hide file tree

Showing 2 changed files with 68 additions and 42 deletions.
diff --git a/oceannavigator/frontend/src/components/AreaWindow.jsx b/oceannavigator/frontend/src/components/AreaWindow.jsx
@@ -407,7 +407,7 @@ class AreaWindow extends React.Component {
         plot_query.contour = this.state.contour;
         plot_query.showarea = this.state.showarea;
         plot_query.variable = this.props.dataset_0.variable;
-        plot_query.projection = this.props.projection;
+        plot_query.projection = this.props.mapSettings.projection;
         plot_query.size = this.state.size;
         plot_query.dpi = this.state.dpi;
         plot_query.interp = this.props.mapSettings.interpType;

diff --git a/plotting/map.py b/plotting/map.py
@@ -136,47 +136,69 @@ def load_data(self):
         width_scale = 1.25
         height_scale = 1.25
 
-        if self.projection == "EPSG:32661":  # north pole projection
-            near_pole, covers_pole = self.pole_proximity(self.points[0])
+        near_pole, covers_pole = self.pole_proximity(self.points[0])
+
+        if self.projection == "EPSG:32661" and (
+            self.centroid[0] > 80 or near_pole or covers_pole
+        ):  # north pole projection
+
             blat = min(self.bounds[0], self.bounds[2])
             blat = 5 * np.floor(blat / 5)
 
-            if self.centroid[0] > 80 or near_pole or covers_pole:
-                self.plot_projection = ccrs.Stereographic(
-                    central_latitude=self.centroid[0],
-                    central_longitude=self.centroid[1],
-                )
-                width_scale = 1.5
-            else:
-                self.plot_projection = ccrs.LambertConformal(
-                    central_latitude=self.centroid[0],
-                    central_longitude=self.centroid[1],
-                )
-        elif self.projection == "EPSG:3031":  # south pole projection
-            near_pole, covers_pole = self.pole_proximity(self.points[0])
+            self.plot_projection = ccrs.Stereographic(
+                central_latitude=90,
+                central_longitude=0,
+            )
+            width_scale = 1.5
+
+            proj_pts = self.plot_projection.transform_points(
+                self.pc_projection,
+                self.points[0][:, 1],
+                self.points[0][:, 0],
+            )
+
+            proj_bounds = np.array(
+                [
+                    [proj_pts[:, 0].min(), proj_pts[:, 1].min()],
+                    [proj_pts[:, 0].max(), proj_pts[:, 1].max()],
+                ]
+            )
+
+        elif self.projection == "EPSG:3031" and (
+            (
+                (self.centroid[0] < -80 or self.bounds[1] < -80 or self.bounds[3] < -80)
+                or covers_pole
+            )
+            or near_pole
+        ):  # south pole projection
             blat = max(self.bounds[0], self.bounds[2])
             blat = 5 * np.ceil(blat / 5)
             # is centerered close to the south pole
-            if (
-                (self.centroid[0] < -80 or self.bounds[1] < -80 or self.bounds[3] < -80)
-                or covers_pole
-            ) or near_pole:
-                self.plot_projection = ccrs.Stereographic(
-                    central_latitude=self.centroid[0],
-                    central_longitude=self.centroid[1],
-                )
-                width_scale = 1.5
-            else:
-                self.plot_projection = ccrs.LambertConformal(
-                    central_latitude=self.centroid[0],
-                    central_longitude=self.centroid[1],
-                )
+
+            self.plot_projection = ccrs.Stereographic(
+                central_latitude=-90,
+                central_longitude=0,
+            )
+            width_scale = 1.5
+
+            proj_pts = self.plot_projection.transform_points(
+                self.pc_projection,
+                self.points[0][:, 1],
+                self.points[0][:, 0],
+            )
+
+            proj_bounds = np.array(
+                [
+                    [proj_pts[:, 0].min(), proj_pts[:, 1].min()],
+                    [proj_pts[:, 0].max(), proj_pts[:, 1].max()],
+                ]
+            )
 
         elif abs(self.centroid[1] - self.bounds[1]) > 90:
 
             if abs(self.bounds[3] - self.bounds[1]) > 360:
                 raise ClientError(
-                    ( # gettext(
+                    (  # gettext(
                         "You have requested an area that exceeds the width \
                         of the world. Thinking big is good but plots need to \
                         be less than 360 deg wide."
@@ -185,21 +207,25 @@ def load_data(self):
             self.plot_projection = ccrs.Mercator(central_longitude=self.centroid[1])
 
         else:
-            self.plot_projection = ccrs.LambertConformal(
-                central_latitude=self.centroid[0], central_longitude=self.centroid[1]
+            self.plot_projection = ccrs.Mercator(
+                central_longitude=self.centroid[1],
+                min_latitude=-85,
+                max_latitude=85,
             )
 
-        proj_bounds = self.plot_projection.transform_points(
-            self.pc_projection,
-            np.array([self.bounds[1], self.bounds[3]]),
-            np.array([self.bounds[0], self.bounds[2]]),
-        )
-        proj_size = np.diff(proj_bounds, axis=0)
+            proj_bounds = self.plot_projection.transform_points(
+                self.pc_projection,
+                np.array([self.bounds[1], self.bounds[3]]),
+                np.array([self.bounds[0], self.bounds[2]]),
+            )
+
+        proj_size = np.diff(proj_bounds[:2, :2], axis=0)
+        proj_centroid = np.mean(proj_bounds[:2, :2], axis=0)
 
         width = proj_size[0][0] * width_scale
         height = proj_size[0][1] * height_scale
 
-        aspect_ratio = height / width
+        aspect_ratio = abs(height / width)
         if aspect_ratio < 1:
             gridx = 500
             gridy = int(500 * aspect_ratio)
@@ -213,8 +239,8 @@ def load_data(self):
             self.plot_projection,
             gridx,
             gridy,
-            x_extents=(-width / 2, width / 2),
-            y_extents=(-height / 2, height / 2),
+            x_extents=(-width / 2, width / 2) + proj_centroid[0],
+            y_extents=(-height / 2, height / 2) + proj_centroid[1],
         )
 
         latlon_grid = self.pc_projection.transform_points(