ESPG fix for 4326 #35

oceanmesh/edgefx.py

@@ -540,7 +540,7 @@ def bathymetric_gradient_sizing_function(
     grid.values = (2 * np.pi / slope_parameter) * np.abs(dp) / (bs + eps)
 
     # Convert back to degrees from meters (if geographic)
-    if crs == "EPSG:4326":
+    if crs == "EPSG:4326" or crs == 4326:
         grid.values /= meters_per_degree
 
     if max_edge_length is not None:
@@ -858,23 +858,27 @@ def wavelength_sizing_function(
     lon, lat = dem.create_grid()
     tmpz = dem.eval((lon, lat))
 
-    if crs == "EPSG:4326":
+    dx, dy = dem.dx, dem.dy  # for gradient function
+
+    if crs == "EPSG:4326" or crs == 4326:
         mean_latitude = np.mean(dem.bbox[2:])
         meters_per_degree = (
             111132.92
             - 559.82 * np.cos(2 * mean_latitude)
             + 1.175 * np.cos(4 * mean_latitude)
             - 0.0023 * np.cos(6 * mean_latitude)
         )
+        dy *= meters_per_degree
+        dx *= meters_per_degree
 
     grid = Grid(
-        bbox=dem.bbox, dx=dem.dx, dy=dem.dy, extrapolate=True, values=0.0, crs=crs
+        bbox=dem.bbox, dx=dx, dy=dy, extrapolate=True, values=0.0, crs=crs
     )
     tmpz[np.abs(tmpz) < 1] = 1
     grid.values = period * np.sqrt(gravity * np.abs(tmpz)) / wl
 
     # Convert back to degrees from meters (if geographic)
-    if crs == "EPSG:4326":
+    if crs == "EPSG:4326" or crs == 4326:
         grid.values /= meters_per_degree
 
     if min_edgelength is None: