Change the projection to better align with x/y

Hard-code the reference latitude
E3SM-Project · Oct 22, 2023 · 8b9b241 · 8b9b241
1 parent cc2a71a
commit 8b9b241
Show file tree

Hide file tree

Showing 4 changed files with 17 additions and 30 deletions.
diff --git a/polaris/ocean/tasks/isomip_plus/mesh/planar.py b/polaris/ocean/tasks/isomip_plus/mesh/planar.py
@@ -52,7 +52,6 @@ def run(self):
         section = config['isomip_plus']
         lx = section.getfloat('lx')
         ly = section.getfloat('ly')
-        lat0 = section.getfloat('lat0')
         buffer = section.getfloat('buffer')
 
         resolution = self.resolution
@@ -72,7 +71,7 @@ def run(self):
         ds_mesh['yIsomipVertex'] = ds_mesh.yVertex
 
         # add latitude and longitude using a stereographic projection
-        projection, lat_lon_projection = get_projections(lat0)
+        projection, lat_lon_projection = get_projections()
         transformer = pyproj.Transformer.from_proj(projection,
                                                    lat_lon_projection)
 

diff --git a/polaris/ocean/tasks/isomip_plus/mesh/spherical.py b/polaris/ocean/tasks/isomip_plus/mesh/spherical.py
@@ -55,9 +55,8 @@ def build_cell_width_lat_lon(self):
 
         lx = section.getfloat('lx')
         ly = section.getfloat('ly')
-        lat0 = section.getfloat('lat0')
         buffer = section.getfloat('buffer')
-        fc = _make_feature(lat0, lx, ly, buffer)
+        fc = _make_feature(lx, ly, buffer)
         fc.to_geojson('isomip_plus_high_res.geojson')
 
         signed_distance = signed_distance_from_geojson(
@@ -80,12 +79,11 @@ def run(self):
         super().run()
 
         ds = xr.open_dataset('base_mesh_without_xy.nc')
-        lat0 = self.config.getfloat('isomip_plus', 'lat0')
-        add_isomip_plus_xy(ds, lat0)
+        add_isomip_plus_xy(ds)
         ds.to_netcdf('base_mesh.nc')
 
 
-def add_isomip_plus_xy(ds, lat0):
+def add_isomip_plus_xy(ds):
     """
     Add x and y coordinates from a stereographic projection to a mesh on a
     sphere
@@ -94,11 +92,8 @@ def add_isomip_plus_xy(ds, lat0):
     ----------
     ds : xarray.Dataset
         The MPAS mesh on a sphere
-
-    lat0 : float
-        The latitude of the origin of the local stereographic project
     """
-    projection, lat_lon_projection = get_projections(lat0)
+    projection, lat_lon_projection = get_projections()
     transformer = pyproj.Transformer.from_proj(lat_lon_projection,
                                                projection)
     lon = np.rad2deg(ds.lonCell.values)
@@ -118,11 +113,11 @@ def add_isomip_plus_xy(ds, lat0):
     ds['yIsomipVertex'] = ('nVertices', y)
 
 
-def _make_feature(lat0, lx, ly, buffer):
+def _make_feature(lx, ly, buffer):
     # a box with a buffer of 80 km surrounding the are of interest
     # (0 <= x <= 800) and (0 <= y <= 80)
     bounds = 1e3 * np.array((-buffer, lx + buffer, -buffer, ly + buffer))
-    projection, lat_lon_projection = get_projections(lat0)
+    projection, lat_lon_projection = get_projections()
     transformer = pyproj.Transformer.from_proj(projection,
                                                lat_lon_projection)
 

diff --git a/polaris/ocean/tasks/isomip_plus/projection.py b/polaris/ocean/tasks/isomip_plus/projection.py
@@ -1,36 +1,30 @@
 import pyproj
 
 
-def get_projection_string(lat0):
+def get_projection_string():
     """
     Get a stereographic projection string that should be used to define
     latitudes and longitudes for the ISOMIP+ x-y coordinates
 
-    Parameters
-    ----------
-    lat0 : float
-        The reference latitude at y=0
-
     Returns
     -------
     proj_str : str
         The projection string
     """
-    proj_str = (f'+proj=stere +lon_0=0 +lat_0={lat0} +lat_ts={lat0} +x_0=0.0 '
-                f'+y_0=0.0 +ellps=WGS84')
+    # y_0 is half way through the domain so the longitude is approximately
+    # like -y
+    #
+    # x_0 is close to lat=-75 and lon=90 so latitude is approximately like x
+    proj_str = ('+proj=stere +lon_0=0 +lat_0=-90 +lat_ts=-75.0 +x_0=-2000e3 '
+                '+y_0=40.0e3 +ellps=WGS84')
     return proj_str
 
 
-def get_projections(lat0):
+def get_projections():
     """
     Get a stereographic and lat-lon projection that can be used to transform
     between lat-lon and ISOMIP+ x-y coordinates
 
-    Parameters
-    ----------
-    lat0 : float
-        The reference latitude at y=0
-
     Returns
     -------
     projection : pyproj.Proj
@@ -39,7 +33,7 @@ def get_projections(lat0):
     lat_lon_projection : pyproj.Proj
         The lat-lon projection
     """
-    projection = pyproj.Proj(get_projection_string(lat0))
+    projection = pyproj.Proj(get_projection_string())
     lat_lon_projection = pyproj.Proj(proj='latlong', datum='WGS84')
 
     return projection, lat_lon_projection
diff --git a/polaris/ocean/tasks/isomip_plus/topo/map.py b/polaris/ocean/tasks/isomip_plus/topo/map.py
@@ -57,9 +57,8 @@ def runtime_setup(self):
         """
         config = self.config
         section = config['isomip_plus']
-        lat0 = section.getfloat('lat0')
         method = section.get('topo_remap_method')
-        proj_str = get_projection_string(lat0)
+        proj_str = get_projection_string()
         self.src_from_proj(filename='input_topo.nc',
                            mesh_name='ISOMIP+_input_topo',
                            x_var='x',