Ignore x0 y0 in dfs2

mikeio/dfs.py

@@ -2,9 +2,10 @@
 from abc import abstractmethod
 from datetime import datetime, timedelta
 from typing import Iterable, List, Optional, Tuple, Union
-
 import numpy as np
 import pandas as pd
+from tqdm import tqdm, trange
+
 from copy import deepcopy
 from mikecore.DfsFactory import DfsFactory
 from mikecore.DfsFile import (
@@ -16,7 +17,7 @@
 )
 from mikecore.DfsFileFactory import DfsFileFactory
 from mikecore.eum import eumQuantity
-from tqdm import tqdm, trange
+from mikecore.Projections import Cartography
 
 from .base import TimeSeries
 from .dataset import Dataset
@@ -751,3 +752,25 @@ def orientation(self):
     def dx(self):
         """Step size in x direction"""
         pass
+
+    def _validate_no_orientation_in_geo(self):
+        if self.is_geo and abs(self._orientation) > 1e-6:
+            raise ValueError("Orientation is not supported for LONG/LAT coordinates")
+
+    def _origin_and_orientation_in_CRS(self):
+        """Project origin and orientation to projected CRS (if not LONG/LAT)"""
+        if self.is_geo:
+            origin = self._longitude, self._latitude
+            orientation = 0.0
+        else:
+            lon, lat = self._longitude, self._latitude
+            cart = Cartography.CreateGeoOrigin(
+                projectionString=self._projstr,
+                lonOrigin=lon,
+                latOrigin=lat,
+                orientation=self._orientation,
+            )
+            # convert origin and orientation to projected CRS
+            origin = tuple(np.round(cart.Geo2Proj(lon, lat), 6))
+            orientation = cart.OrientationProj
+        return origin, orientation

mikeio/dfs2.py

@@ -1,14 +1,14 @@
 import os
 from copy import deepcopy
-
 import numpy as np
 import pandas as pd
+from tqdm import tqdm
+
 from mikecore.DfsFactory import DfsBuilder, DfsFactory
 from mikecore.DfsFile import DfsFile, DfsSimpleType
 from mikecore.DfsFileFactory import DfsFileFactory
 from mikecore.eum import eumQuantity, eumUnit
 from mikecore.Projections import Cartography
-from tqdm import tqdm
 
 from . import __dfs_version__
 from .dataset import Dataset
@@ -44,14 +44,14 @@ def _write_dfs2_header(filename, ds: Dataset, title="") -> DfsFile:
 
     factory = DfsFactory()
     _write_dfs2_spatial_axis(builder, factory, geometry)
-    proj = geometry.projection_string
+    proj_str = geometry.projection_string
     origin = geometry.origin
     orient = geometry.orientation
 
     if geometry.is_geo:
-        proj = factory.CreateProjectionGeoOrigin(proj, *origin, orient)
+        proj = factory.CreateProjectionGeoOrigin(proj_str, *origin, orient)
     else:
-        cart: Cartography = Cartography.CreateProjOrigin(proj, *origin, orient)
+        cart: Cartography = Cartography.CreateProjOrigin(proj_str, *origin, orient)
         proj = factory.CreateProjectionGeoOrigin(
             wktProjectionString=geometry.projection,
             lon0=cart.LonOrigin,
@@ -105,62 +105,35 @@ class Dfs2(_Dfs123):
 
     _ndim = 2
 
-    def __init__(self, filename=None, type="horizontal"):
+    def __init__(self, filename=None, type: str = "horizontal"):
         super().__init__(filename)
 
         self._dx = None
         self._dy = None
         self._nx = None
         self._ny = None
-        self._x0 = 0
-        self._y0 = 0
+        self._x0 = 0.0
+        self._y0 = 0.0
         self.geometry = None
 
         if filename:
-            self._read_dfs2_header()
-            is_spectral = type == "spectral"
-
-            if self._projstr == "LONG/LAT":
-                if np.abs(self._orientation) < 1e-6:
-                    origin = self._longitude, self._latitude
-                    self.geometry = Grid2D(
-                        dx=self._dx,
-                        dy=self._dy,
-                        nx=self._nx,
-                        ny=self._ny,
-                        x0=self._x0,
-                        y0=self._y0,
-                        origin=origin,
-                        projection=self._projstr,
-                        is_spectral=is_spectral,
-                    )
-                else:
-                    raise ValueError(
-                        "LONG/LAT with non-zero orientation is not supported"
-                    )
-            else:
-                lon, lat = self._longitude, self._latitude
-                cart = Cartography.CreateGeoOrigin(
-                    projectionString=self._projstr,
-                    lonOrigin=self._longitude,
-                    latOrigin=self._latitude,
-                    orientation=self._orientation,
-                )
-                origin_projected = np.round(cart.Geo2Proj(lon, lat), 7)
-                orientation_projected = cart.OrientationProj
-
-                self.geometry = Grid2D(
-                    dx=self._dx,
-                    dy=self._dy,
-                    nx=self._nx,
-                    ny=self._ny,
-                    x0=self._x0,
-                    y0=self._y0,
-                    orientation=orientation_projected,
-                    origin=origin_projected,
-                    projection=self._projstr,
-                    is_spectral=is_spectral,
-                )
+            is_spectral = type.lower() in ["spectral", "spectra", "spectrum"]
+            self._read_dfs2_header(read_x0y0=is_spectral)
+            self._validate_no_orientation_in_geo()
+            origin, orientation = self._origin_and_orientation_in_CRS()
+
+            self.geometry = Grid2D(
+                dx=self._dx,
+                dy=self._dy,
+                nx=self._nx,
+                ny=self._ny,
+                x0=self._x0,
+                y0=self._y0,
+                orientation=orientation,
+                origin=origin,
+                projection=self._projstr,
+                is_spectral=is_spectral,
+            )
 
     def __repr__(self):
         out = ["<mikeio.Dfs2>"]
@@ -185,14 +158,15 @@ def __repr__(self):
 
         return str.join("\n", out)
 
-    def _read_dfs2_header(self):
+    def _read_dfs2_header(self, read_x0y0: bool = False):
         if not os.path.isfile(self._filename):
             raise Exception(f"file {self._filename} does not exist!")
 
         self._dfs = DfsFileFactory.Dfs2FileOpen(self._filename)
         self._source = self._dfs
-        self._x0 = self._dfs.SpatialAxis.X0
-        self._y0 = self._dfs.SpatialAxis.Y0
+        if read_x0y0:
+            self._x0 = self._dfs.SpatialAxis.X0
+            self._y0 = self._dfs.SpatialAxis.Y0
         self._dx = self._dfs.SpatialAxis.Dx
         self._dy = self._dfs.SpatialAxis.Dy
         self._nx = self._dfs.SpatialAxis.XCount

mikeio/dfs3.py

@@ -2,11 +2,13 @@
 
 import numpy as np
 import pandas as pd
+
 from mikecore.DfsBuilder import DfsBuilder
 from mikecore.DfsFactory import DfsFactory
 from mikecore.DfsFile import DfsFile, DfsSimpleType
 from mikecore.DfsFileFactory import DfsFileFactory
 from mikecore.eum import eumQuantity, eumUnit
+from mikecore.Projections import Cartography
 
 from . import __dfs_version__
 from .dataset import Dataset
@@ -34,23 +36,23 @@ def _write_dfs3_header(filename, ds: Dataset, title="") -> DfsFile:
 
     factory = DfsFactory()
     _write_dfs3_spatial_axis(builder, factory, geometry)
-    origin = geometry._origin  # Origin in geographical coordinates
-    orient = geometry._orientation
+    origin = geometry.origin  # Origin in geographical coordinates
+    orient = geometry.orientation
 
-    # if geometry.is_geo:
-    proj = factory.CreateProjectionGeoOrigin(
-        geometry.projection_string, *origin, orient
-    )
-    # else:
-    #    cart: Cartography = Cartography.CreateProjOrigin(
-    #        geometry.projection_string, *origin, orient
-    #    )
-    #    proj = factory.CreateProjectionGeoOrigin(
-    #        wktProjectionString=geometry.projection,
-    #        lon0=cart.LonOrigin,
-    #        lat0=cart.LatOrigin,
-    #        orientation=cart.Orientation,
-    #    )
+    if geometry.is_geo:
+        proj = factory.CreateProjectionGeoOrigin(
+            geometry.projection_string, *origin, orient
+        )
+    else:
+        cart: Cartography = Cartography.CreateProjOrigin(
+            geometry.projection_string, *origin, orient
+        )
+        proj = factory.CreateProjectionGeoOrigin(
+            wktProjectionString=geometry.projection,
+            lon0=cart.LonOrigin,
+            lat0=cart.LatOrigin,
+            orientation=cart.Orientation,
+        )
 
     builder.SetGeographicalProjection(proj)
 
@@ -110,13 +112,16 @@ def __init__(self, filename=None):
         self._nx = None
         self._ny = None
         self._nz = None
-        self._x0 = 0
-        self._y0 = 0
-        self._z0 = 0
+        self._x0 = 0.0
+        self._y0 = 0.0
+        self._z0 = 0.0
         self.geometry = None
 
         if filename:
             self._read_dfs3_header()
+            self._validate_no_orientation_in_geo()
+            origin, orientation = self._origin_and_orientation_in_CRS()
+
             self.geometry = Grid3D(
                 x0=self._x0,
                 dx=self._dx,
@@ -127,9 +132,9 @@ def __init__(self, filename=None):
                 z0=self._z0,
                 dz=self._dz,
                 nz=self._nz,
-                origin=(self._longitude, self._latitude),
+                origin=origin,
                 projection=self._projstr,
-                orientation=self._orientation,
+                orientation=orientation,
             )
 
     def __repr__(self):
@@ -154,17 +159,19 @@ def __repr__(self):
 
         return str.join("\n", out)
 
-    def _read_dfs3_header(self):
+    def _read_dfs3_header(self, read_x0y0z0: bool = False):
         if not os.path.isfile(self._filename):
             raise Exception(f"file {self._filename} does not exist!")
 
         self._dfs = DfsFileFactory.Dfs3FileOpen(self._filename)
 
-        self.source = self._dfs
+        self._source = self._dfs
+
+        if read_x0y0z0:
+            self._x0 = self._dfs.SpatialAxis.X0
+            self._y0 = self._dfs.SpatialAxis.Y0
+            self._z0 = self._dfs.SpatialAxis.Z0
 
-        self._x0 = self._dfs.SpatialAxis.X0
-        self._y0 = self._dfs.SpatialAxis.Y0
-        self._z0 = self._dfs.SpatialAxis.Z0
         self._dx = self._dfs.SpatialAxis.Dx
         self._dy = self._dfs.SpatialAxis.Dy
         self._dz = self._dfs.SpatialAxis.Dz
@@ -406,3 +413,8 @@ def dz(self):
     @property
     def shape(self):
         return (self._n_timesteps, self._nz, self._ny, self._nx)
+
+    @property
+    def is_geo(self):
+        """Are coordinates geographical (LONG/LAT)?"""
+        return self._projstr == "LONG/LAT"