inverse distance weighting

GEDI_11SLT.jpeg.aux.xml

@@ -0,0 +1,24 @@
+<PAMDataset>
+  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]</SRS>
+  <GeoTransform>  3.0000000000000000e+05,  3.0000000000000000e+01,  0.0000000000000000e+00,  3.8000400000000000e+06,  0.0000000000000000e+00, -3.0000000000000000e+01</GeoTransform>
+  <Metadata domain="IMAGE_STRUCTURE">
+    <MDI key="SOURCE_COLOR_SPACE">YCbCr</MDI>
+    <MDI key="INTERLEAVE">PIXEL</MDI>
+    <MDI key="COMPRESSION">JPEG</MDI>
+  </Metadata>
+  <PAMRasterBand band="1">
+    <Metadata domain="IMAGE_STRUCTURE">
+      <MDI key="COMPRESSION">JPEG</MDI>
+    </Metadata>
+  </PAMRasterBand>
+  <PAMRasterBand band="2">
+    <Metadata domain="IMAGE_STRUCTURE">
+      <MDI key="COMPRESSION">JPEG</MDI>
+    </Metadata>
+  </PAMRasterBand>
+  <PAMRasterBand band="3">
+    <Metadata domain="IMAGE_STRUCTURE">
+      <MDI key="COMPRESSION">JPEG</MDI>
+    </Metadata>
+  </PAMRasterBand>
+</PAMDataset>

rasters/local_UTM_proj4_from_lat_lon.py

@@ -15,6 +15,6 @@ def local_UTM_proj4_from_lat_lon(lat: float, lon: float) -> str:
         str: The proj4 string for the local UTM projection.
     """
     UTM_zone = (np.floor((lon + 180) / 6) % 60) + 1
-    UTM_proj4 = f"+proj=utm +zone={UTM_zone} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs"
+    UTM_proj4 = f"+proj=utm +zone={int(UTM_zone)} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs"
 
     return UTM_proj4

rasters/point.py

@@ -1,15 +1,18 @@
 from __future__ import annotations
 
-from typing import Union
+from typing import Union, Iterable
 
 import geopandas as gpd
 import shapely
 from shapely.geometry.base import CAP_STYLE, JOIN_STYLE
 
+from rasters.wrap_geometry import wrap_geometry
+
 from .constants import *
 from .CRS import CRS, WGS84
 from .vector_geometry import VectorGeometry, SingleVectorGeometry
 from .polygon import Polygon
+from .bbox import BBox
 
 class Point(SingleVectorGeometry):
     """
@@ -128,3 +131,69 @@ def buffer(
         )
         # Return a new Polygon object with the buffer geometry and the point's CRS
         return Polygon(buffer_geom, crs=self.crs)
+
+    @property
+    def bbox(self) -> BBox:
+        """
+        Returns the bounding box of the point.
+
+        Returns:
+            BBox: The bounding box of the point.
+        """
+        return BBox(self.x, self.y, self.x, self.y, crs=self.crs)
+
+    def distance(self, other: Point) -> float:
+        """
+        Returns the distance between this point and another point.
+
+        Args:
+            other (Point): The other point.
+
+        Returns:
+            float: The distance between the two points.
+        """
+        if isinstance(other, shapely.geometry.Point):
+            other = Point(other, crs=WGS84)
+
+        projected = self.projected
+        # print(f"projects CRS: {projected.crs}")
+        # print(f"other: {other} {type(other)} {other.crs}")
+        other_projected = other.to_crs(projected.crs)
+        other_projected._crs = projected.crs
+        # print(f"other_projected.crs: {other_projected.crs}")
+        # print(f"other_projected: {other_projected} {type(other_projected)}")
+        x, y = projected.geometry.x, projected.geometry.y
+        # print(f"x: {x}, y: {y}")
+        other_x, other_y = other.geometry.x, other.geometry.y
+        # print(f"other_x: {other_x}, other_y: {other_y}")
+        distance = ((x - other_x) ** 2 + (y - other_y) ** 2) ** 0.5
+
+        return distance
+
+    from .multi_point import MultiPoint
+
+    def distances(self, points: Union[MultiPoint, Iterable[Point]]) -> gpd.GeoDataFrame:
+        """
+        Calculate the distances between this geometry and a set of points.
+
+        Args:
+            points (Union[MultiPoint, Iterable[Point]]): The points to calculate distances to.
+
+        Returns:
+            gpd.GeoDataFrame: A GeoDataFrame with the distances between the geometry and the points.
+                                The GeoDataFrame has a 'distance' column and a 'geometry' column
+                                containing LineStrings between the geometry and each point.
+        """
+        points = wrap_geometry(points)
+
+        geometry = []
+        distance_column = []
+
+        for point in [wrap_geometry(point) for point in points.geoms]:
+            point._crs = points._crs
+            geometry.append(shapely.geometry.LineString([self.geometry, point.geometry]))
+            distance_column.append(self.distance(point))
+
+        distances = gpd.GeoDataFrame({"distance": distance_column}, geometry=geometry)
+
+        return distances

rasters/raster.py

@@ -1147,13 +1147,18 @@ def pixel_outlines(self) -> gpd.GeoDataFrame:
         return gpd.GeoDataFrame({"value": self.array.flatten()}, geometry=list(self.geometry.pixel_outlines))
 
     def IDW(self, geometry: VectorGeometry, power: float = 2):
-        if isinstance(geometry, shapely.geometry.Point):
+        from .point import Point
+        from .vector_geometry import SingleVectorGeometry
+
+        if isinstance(geometry, (Point, shapely.geometry.Point)):
             geometry = wrap_geometry(geometry)
             pixel_centroids = self.geometry.pixel_centroids
             distances = geometry.distances(pixel_centroids)
             value = self.array.flatten()
             distance = np.array(distances["distance"])
+            # print(f"distances: {distance}")
             weight = 1 / distance ** power
+            # print(f"weights: {weight}")
             weighted = value * weight
             result = np.nansum(weighted) / np.nansum(weight)
 
@@ -1278,7 +1283,7 @@ def to_geotiff(
             filename: str,
             compression: str = None,
             preview_quality: int = 20,
-            cmap: Union[Colormap, str] = DEFAULT_CMAP,
+            cmap: Union[Colormap, str] = None,
             use_compression: bool = True,
             include_preview: bool = True,
             overwrite: bool = True,

rasters/raster_geometry.py

@@ -1087,6 +1087,7 @@ def pixel_centroids(self) -> MultiPoint:
 
         x, y = self.xy
         pixel_centroids = wrap_geometry(MultiPoint(np.stack([x.flatten(), y.flatten()], axis=1)), crs=self.crs)
+        pixel_centroids._crs = self.crs
 
         return pixel_centroids
 

rasters/vector_geometry.py

@@ -145,31 +145,6 @@ def projected(self) -> VectorGeometry:
         else:
             return self
 
-    def distances(self, points: Union[MultiPoint, Iterable[Point]]) -> gpd.GeoDataFrame:
-        """
-        Calculate the distances between this geometry and a set of points.
-
-        Args:
-            points (Union[MultiPoint, Iterable[Point]]): The points to calculate distances to.
-
-        Returns:
-            gpd.GeoDataFrame: A GeoDataFrame with the distances between the geometry and the points.
-                                The GeoDataFrame has a 'distance' column and a 'geometry' column
-                                containing LineStrings between the geometry and each point.
-        """
-        points = [wrap_geometry(point) for point in points]
-
-        geometry = []
-        distance_column = []
-
-        for point in points:
-            geometry.append(shapely.geometry.LineString([self, point]))
-            distance_column.append(self.projected.distance(point.projected))
-
-        distances = gpd.GeoDataFrame({"distance": distance_column}, geometry=geometry)
-
-        return distances
-
 
 class SingleVectorGeometry(VectorGeometry):
     """
@@ -181,4 +156,19 @@ class MultiVectorGeometry(VectorGeometry):
     """
     Base class for multi vector geometries (e.g., MultiPoint, MultiLineString, MultiPolygon).
     """
-    pass
+    @property
+    def geoms(self):
+        return self.geometry.geoms
+
+    def __getitem__(self, index):
+        if isinstance(index, int):
+            return self.geoms[index]
+
+        return self.geometry[index]
+
+    # @property
+    # def geoms(self) -> list[VectorGeometry]:
+    #     """
+    #     The list of geometries in the multi geometry.
+    #     """ 
+    #     return [self.contain(geometry) for geometry in self.geometry.geoms]

rasters/wrap_geometry.py

@@ -36,6 +36,10 @@ def wrap_geometry(geometry: Any, crs: Union[CRS, str] = None) -> SpatialGeometry
     Raises:
         ValueError: If the geometry type is not supported.
     """
+    from .point import Point
+    from .multi_point import MultiPoint
+    from .polygon import Polygon
+    from .multi_polygon import MultiPolygon
 
     if isinstance(geometry, SpatialGeometry):
         # If the geometry is already a SpatialGeometry, return it as is