Merge pull request #50 from ruithnadsteud/crs

Polygon and dataset CRS objects kwarg addition

.circleci/config.yml

@@ -199,4 +199,4 @@ workflows:
        - tests:
            name: "Python 3.8 tests"
            tag: "3.8.7"
-           ytdev: 1
+           ytdev: 1

setup.py

@@ -81,6 +81,7 @@ def run(self):
     packages=["yt_georaster"],
     keywords=["GeoTiff", "GTiff", "raster"],
     install_requires=[
+        "cython",
         "fiona",
         "gdal",
         "numpy",

yt_georaster/data_structures.py

@@ -3,6 +3,7 @@
 import rasterio
 from rasterio import warp
 from rasterio.windows import from_bounds
+from rasterio.crs import CRS
 import re
 import weakref
 
@@ -56,8 +57,8 @@ def __init__(self, gridobj, left_edge, right_edge):
         self._parent_id = -1
         self.Level = 0
 
-        self.LeftEdge = self.ds.arr(left_edge, "m")
-        self.RightEdge = self.ds.arr(right_edge, "m")
+        self.LeftEdge = self.ds.arr(left_edge, self.ds.parameters['units'])
+        self.RightEdge = self.ds.arr(right_edge, self.ds.parameters['units'])
         # Make sure z dimension edges are the same as parent grid.
         self.LeftEdge[2] = gridobj.LeftEdge[2]
         self.RightEdge[2] = gridobj.RightEdge[2]
@@ -292,10 +293,11 @@ class GeoRasterDataset(Dataset):
     refine_by = 2
     _con_attrs = ()
 
-    def __init__(self, *args, field_map=None):
+    def __init__(self, *args, field_map=None, crs=None):
         self.filename_list = args
         filename = args[0]
         self.field_map = field_map
+        self.crs = crs
         super().__init__(filename, self._dataset_type, unit_system="mks")
         self.data = self.index.grids[0]
         self._added_fields = []
@@ -324,9 +326,71 @@ def _parse_parameter_file(self):
                 v = f.meta[key]
                 self.parameters[key] = v
             self.parameters["res"] = f.res
-            self.parameters["src"] = f.crs
         self.current_time = 0
 
+        # overwrite crs if one is provided by user
+        if not (self.crs is None):
+            # make sure user provided CRS is valid CRS object
+            if not isinstance(self.crs, CRS):
+                if isinstance(self.crs, int):
+                    # assume epsg number
+                    self.crs = CRS.from_epsg(self.crs)
+                elif isinstance(self.crs, dict):
+                    self.crs = CRS.from_dict(**self.crs)
+                else:
+                    self.crs = CRS.from_string(self.crs)
+
+            # get reprojected transform
+            left_edge = self.parameters["transform"] * (0, 0)
+            right_edge = self.parameters["transform"] * (
+                self.parameters["width"],
+                self.parameters["height"]
+            )
+            transform, width, height = warp.calculate_default_transform(
+                self.parameters["crs"],
+                self.crs,
+                self.parameters["width"],
+                self.parameters["height"],
+                left=left_edge[0],
+                bottom=left_edge[1],
+                right=right_edge[0],
+                top=right_edge[1],
+                # resolution=self.parameters["transform"][0]
+                dst_width=self.parameters["width"],
+                dst_height=self.parameters["height"]
+            )  # current solution can create rectangular pixels
+            # xs, ys = warp.transform(
+            #     self.parameters["crs"],
+            #     dst_crs,
+            #     [left_edge[0], right_edge[0]],
+            #     [left_edge[1], right_edge[1]],
+            #     zs=None
+            # )
+            # update parameters
+            self.parameters["res"] = (transform[0], -transform[4])
+            self.parameters["crs"] = self.crs
+            self.parameters["transform"] = transform
+            self.parameters["width"] = width
+            self.parameters["height"] = height
+        else:
+            # if no crs has be provided replace None with base image CRS
+            self.crs = self.parameters["crs"]
+
+        # get units and conversion factor to metres
+        self.parameters["units"] = self.parameters["crs"].linear_units
+        # for non-projected crs this is unknown
+        if self.parameters["units"] == 'unknown':
+            mylog.warning(
+                f"Dataset CRS {self.parameters['crs']} "
+                "units are 'unknown'. Using meters."
+            )
+            self.parameters["units"] = 'm'  # just a place holder
+        else:
+            mylog.info(
+                f"Dataset CRS {self.parameters['crs']} "
+                f"units are '{self.parameters['units']}'. "
+            )
+        # set domain
         width = self.parameters["width"]
         height = self.parameters["height"]
         transform = self.parameters["transform"]
@@ -337,9 +401,18 @@ def _parse_parameter_file(self):
         rast_right = np.concatenate([transform * (width, height), [1]])
         # save dimensions that need to be flipped
         self._flip_axes = np.where(rast_left > rast_right)[0]
-        self.domain_left_edge = self.arr(np.min([rast_left, rast_right], axis=0), "m")
-        self.domain_right_edge = self.arr(np.max([rast_left, rast_right], axis=0), "m")
-        self.resolution = self.arr(self.parameters["res"], "m")
+        self.domain_left_edge = self.arr(
+            np.min([rast_left, rast_right], axis=0),
+            self.parameters["units"]
+        )
+        self.domain_right_edge = self.arr(
+            np.max([rast_left, rast_right], axis=0),
+            self.parameters["units"]
+        )
+        self.resolution = self.arr(
+            self.parameters["res"],
+            self.parameters["units"]
+        )
 
     def _setup_classes(self):
         super()._setup_classes()
@@ -649,8 +722,8 @@ def __init__(self, parent_ds, left_edge, right_edge):
         self._parent_ds = parent_ds
         self._index_class = parent_ds._index_class
         self._dataset_type = parent_ds._dataset_type
-        self.domain_left_edge = parent_ds.arr(left_edge, "m")
-        self.domain_right_edge = parent_ds.arr(right_edge, "m")
+        self.domain_left_edge = parent_ds.arr(left_edge, parent_ds.parameters["units"])
+        self.domain_right_edge = parent_ds.arr(right_edge, parent_ds.parameters["units"])
         self.domain_dimensions = (
             parent_ds.domain_dimensions
             * (self.domain_right_edge - self.domain_left_edge)

yt_georaster/image_types.py

@@ -34,7 +34,11 @@ def identify(self, filename):
 
 
 class Sentinel2(SatGeoImage):
-    _regex = re.compile(r"([A-Za-z0-9]+_[A-Za-z0-9]+)_([A-Za-z0-9]+)(?:_\d+m)?$")
+    _regex = re.compile(
+        r'^[A-Za-z0-9]+_[A-Za-z0-9]+_([A-Za-z0-9]+)_[A-Za-z0-9]+_'
+        r'[A-Za-z0-9]+_[A-Za-z0-9]+_[A-Za-z0-9]+_'
+        r'([A-Za-z0-9]+)(?:_\\d+m)?$'
+    )
     _suffix = "jp2"
     _field_prefix = "S2"
     _band_aliases = (

yt_georaster/io.py

@@ -106,11 +106,14 @@ def _read_rasterio_data(self, selector, grid, field):
 
         # Resample to base resolution if necessary.
         image_resolution = src.res[0]
+        image_units = src.crs.linear_units
         base_resolution = self.ds.resolution.d[0]
+        base_units = self.ds.parameters['units']
         if image_resolution != base_resolution:
             scale_factor = image_resolution / base_resolution
             mylog.info(
-                f"Resampling {field}: {image_resolution} to {base_resolution} m."
+                f"Resampling {field}: {image_resolution} {image_units} "
+                f"to {base_resolution} {base_units}."
             )
             data = zoom(data, scale_factor, order=0)
 

yt_georaster/polygon.py

@@ -8,6 +8,8 @@
 import numpy as np
 from shapely.geometry import Polygon, MultiPolygon
 from shapely.ops import unary_union
+from rasterio.crs import CRS
+from rasterio.warp import transform_geom
 
 from yt_georaster.polygon_selector import PolygonSelector
 
@@ -40,32 +42,65 @@ class YTPolygon(YTSelectionContainer3D):
     _con_args = ("filename",)
 
     # add more arguments, like path to a shape file or a shapely Polygon object
-    def __init__(self, filename, ds=None, field_parameters=None):
+    def __init__(self, filename, ds=None, field_parameters=None, crs=None):
         validate_object(ds, Dataset)
         validate_object(field_parameters, dict)
-
-        self.filename = filename
-
-        # read shapefile with fiona
-        with fiona.open(filename, "r") as shapefile:
-            shapes_from_file = [feature["geometry"] for feature in shapefile]
-
-        mylog.info(f"Number of features in file: {len(shapes_from_file)}")
-
-        # save number of polygons
-        self._number_features = len(shapes_from_file)
-
-        # convert all polyogn features in shapefile to list of shapely polygons
-        polygons = [
-            Polygon(shapes_from_file[i]["coordinates"][0])
-            for i in range(len(shapes_from_file))
-        ]
-
-        #  fix invalid MultiPolygons
-        m = MultiPolygon(polygons)
-
-        # join all shapely polygons to a single layer
-        self.polygon = unary_union(m)
+        self.src_crs = crs
+        if isinstance(filename, str):
+            self.filename = filename
+
+            # read shapefile with fiona
+            with fiona.open(filename, "r") as shapefile:
+                shapes_from_file = [feature["geometry"] for feature in shapefile]
+                self.src_crs = CRS.from_dict(**shapefile.crs)  # shapefile crs
+
+            # save number of polygons
+            self._number_features = len(shapes_from_file)
+
+            # reproject to datasets crs
+            for i in range(self._number_features):
+                shapes_from_file[i] = transform_geom(
+                    f'EPSG:{self.src_crs.to_epsg()}',
+                    f'EPSG:{ds.parameters["crs"].to_epsg()}',
+                    shapes_from_file[i]
+                )
+            # convert all polygon features in shapefile to list of shapely polygons
+            polygons = [
+                Polygon(shapes_from_file[i]["coordinates"][0])
+                for i in range(self._number_features)
+            ]
+            # fix invalid MultiPolygons
+            m = MultiPolygon(polygons)
+            # join all shapely polygons to a single layer
+            self.polygon = unary_union(m)
+
+        elif isinstance(filename, Polygon):
+            # only one polygon
+            self._number_features = 1
+            self.polygon = filename
+            if not (self.src_crs is None):
+                self._reproject_polygon(ds.parameters['crs'])
+
+        elif isinstance(filename, MultiPolygon):
+            # only one polygon
+            self._number_features = len(filename.geoms)
+            self.polygon = unary_union(filename)
+            if not (self.src_crs is None):
+                self._reproject_polygon(ds.parameters['crs'])
+
+        elif isinstance(filename, list):
+            # assume list of shapely polygons
+            self._number_features = len(filename)
+            # fix invalid MultiPolygons
+            m = MultiPolygon(filename)
+            # join all shapely polygons to a single layer
+            self.polygon = unary_union(m)
+            if not (self.src_crs is None):
+                self._reproject_polygon(ds.parameters['crs'])
+
+        mylog.info(
+            f"Number of features in poly object: {self._number_features}"
+        )
 
         # define coordinates of center
         self.center = [
@@ -76,6 +111,60 @@ def __init__(self, filename, ds=None, field_parameters=None):
         data_source = None
         super().__init__(self.center, ds, field_parameters, data_source)
 
+    def _reproject_polygon(self, dst_crs):
+        """
+        Reproject polygon objects to destination projection.
+
+        fiona transform geom object does not handle shapely objects but
+        "GeoJSON-like" dictionaries instead.
+        """
+        # get epsg
+        src_crs = self.src_crs
+        if isinstance(src_crs, int):
+            # assume epsg number
+            epsg = src_crs
+            src_crs = CRS.from_epsg(epsg)
+        elif isinstance(src_crs, dict):
+            src_crs = CRS.from_dict(**src_crs)
+            epsg = src_crs.to_epsg()
+        elif isinstance(src_crs, CRS):
+            epsg = src_crs.to_epsg()
+        else:
+            src_crs = CRS.from_string(src_crs)
+            epsg = src_crs.to_epsg()
+
+        if dst_crs.to_epsg() != epsg:
+            if isinstance(self.polygon, MultiPolygon):
+                coords = [list(geom.exterior.coords) for geom in self.polygon.geoms]
+                type_str = 'MultiPolygon'
+                transformed_poly = transform_geom(
+                    src_crs,
+                    dst_crs,
+                    {
+                        'type': type_str,
+                        'coordinates': coords
+                    }
+                )
+                # convert all polygon features back to shapely polygons
+                polygons = [
+                    Polygon(transformed_poly[i]["coordinates"][0])
+                    for i in range(len(transformed_poly))
+                ]
+                # join all shapely polygons to a single layer
+                self.polygon = unary_union(MultiPolygon(polygons))
+            else:
+                coords = [list(self.polygon.exterior.coords)]
+                type_str = 'Polygon'
+                transformed_poly = transform_geom(
+                    src_crs,
+                    dst_crs,
+                    {
+                        'type': type_str,
+                        'coordinates': coords
+                    }
+                )
+                self.polygon = Polygon(transformed_poly['coordinates'][0])
+
     def _get_bbox(self):
         """
         Return the minimum bounding box for the polygon.