Merge pull request #17 from Open-Source-Agriculture/GeoTiffObject

Bounding box features

README.md

@@ -18,7 +18,7 @@ pip install geotiff
 
 ### Usage
 
-Read the GeoTiff to an array
+#### Read the GeoTiff to an array
 
 ```python
 from geotiff import GeoTiff
@@ -27,28 +27,55 @@ geoTiff = GeoTiff(tiff_file)
 array = geoTiff.read()
 ```
 
-Read a sections of a large tiff using a WGS 84 area
+Get bounding box info about the tiff
+
+```python
+# in the original CRS
+geotiff.tif_bBox
+# as WGS 84
+geotiff.tif_bBox_wgs_84
+```
+
+Get coordinates of a point/pixel
+
+```python
+i=5
+j=6
+geoTiff.get_wgs_84_coords(i, j)
+```
+
+Get the original crs code
+
+```python
+geotiff.crs_code
+```
+
+#### Read a sections of a large tiff using a WGS 84 area
 
 ```python
 from geotiff import GeoTiff
 
+# in WGS 84
 area_box = [(138.632071411, -32.447310785), (138.644218874, -32.456979174)]
 geotiff = GeoTiff(tiff_file)
 array = geotiff.read_box(area_box)
 ```
 
-This will detect and convert coordinates into WGS 84
-
-You can also get the coordinated of the actual bounding box
+Getting bounding box information
 
 ```python
-geotiff.tif_bBox
+# col and row indexes of the cut area
+int_box = geoTiff.get_int_box(area_box)
+# lon and lat coords of the cut points/pixels
+geoTiff.get_bBox_wgs_84(area_box)
 ```
 
-And the original crs code
+Get coordinates of a point/pixel
 
 ```python
-geotiff.crs_code
+i=int_box[0][0] + 5
+j=int_box[0][1] + 6
+geoTiff.get_wgs_84_coords(i, j)
 ```
 
 ### Contributing

example.py

@@ -4,19 +4,25 @@
 from geotiff import GeoTiff # type: ignore
 
 
+filename = "dem.tif"
 filename = "dem.tif"
 dir = "./tests/inputs/"
 tiff_file = os.path.join(dir, filename)
+# tiff_file = "/home/kipling/Documents/fsm_sample_data/ntz52.tif"
 # 138.632071411 -32.447310785 138.644218874 -32.456979174
 bounding_box: List[Tuple[float, float]] = [(138.632071411, -32.447310785), (138.644218874, -32.456979174)]
 
 
 print("testing read tiff")
 print(f"reading: {tiff_file}")
 print(f"Using bBox: {bounding_box}")
-array = GeoTiff(tiff_file).read_box(bounding_box)
+geotiff: GeoTiff = GeoTiff(tiff_file)
+print(geotiff.tif_bBox)
+print(geotiff.crs_code)
+array = geotiff.read_box(bounding_box)
 print("Sample array:")
 print(array)
 print(array.shape)
 assert isinstance(array, np.ndarray)
+print(geotiff.tif_bBox_wgs_84)
 print("end")

geotiff/__init__.py

@@ -155,8 +155,7 @@ def _convert_to_wgs_84(self, crs_code: int, xxyy: Tuple[float, float]) -> Tuple[
         yy: float = xxyy[1]
         crs_4326: CRS = CRS("WGS84")
         crs_proj: CRS = CRS.from_epsg(crs_code)
-        transformer: Transformer = Transformer.from_crs(
-            crs_proj, crs_4326, always_xy=True)
+        transformer: Transformer = Transformer.from_crs(crs_proj, crs_4326, always_xy=True)
         return(transformer.transform(xx, yy))
 
 
@@ -169,6 +168,7 @@ def _convert_from_wgs_84(self, crs_code: int, xxyy: Tuple[float, float]) -> Tupl
             crs_4326, crs_proj, always_xy=True)
         return(transformer.transform(xx, yy))
 
+
     def _get_x_int(self, lon) -> int:
         step_x: float = float(
             self.tifShape[1]/(self.tif_bBox[1][0] - self.tif_bBox[0][0]))
@@ -178,6 +178,28 @@ def _get_y_int(self, lat) -> int:
         step_y: float = self.tifShape[0]/(self.tif_bBox[1][1] - self.tif_bBox[0][1])
         return(int(step_y*(lat - self.tif_bBox[0][1])))
 
+    def get_wgs_84_coords(self, i: int, j)-> Tuple[float, float]:
+        """for a given i, j in the entire tiff array,
+        returns the wgs_84 coordinates
+
+        Args:
+            i (int): col number of the array
+            j (int): row number of the array
+
+        Returns:
+            Tuple[float, float]: lon, lat
+        """
+        x, y = self.tifTrans.get_xy(i, j)
+        return(self._convert_to_wgs_84(self.crs_code, (x, y)))
+
+
+    @property
+    def tif_bBox_wgs_84(self) -> BBox: 
+        right_top = self._convert_to_wgs_84(self.crs_code, self.tif_bBox[0])
+        left_bottom = self._convert_to_wgs_84(self.crs_code, self.tif_bBox[1])
+        return(right_top, left_bottom)
+
+
     def get_int_box(self, bBox: BBox) -> BBoxInt:
         """Gets the intiger array index values based on a bounding box
 
@@ -192,16 +214,16 @@ def get_int_box(self, bBox: BBox) -> BBoxInt:
             BBoxInt: array index values
         """
         b_bBox = (self._convert_from_wgs_84(self.crs_code, bBox[0]), self._convert_from_wgs_84(self.crs_code, bBox[1]))
-        x_min: int = self._get_x_int(b_bBox[0][0])
-        y_min: int = self._get_y_int(b_bBox[0][1])
+        x_min: int = self._get_x_int(b_bBox[0][0]) + int(not self.tif_bBox[0][0]==b_bBox[0][0])
+        y_min: int = self._get_y_int(b_bBox[0][1]) + int(not self.tif_bBox[0][1]==b_bBox[0][1])
         x_max: int = self._get_x_int(b_bBox[1][0])
         y_max: int = self._get_y_int(b_bBox[1][1])
 
-        shp_bBox = [self.tifTrans.get_xy(x_min,y_min),  self.tifTrans.get_xy(x_max+1,y_max+1)]
-        check = (shp_bBox[0][0] < b_bBox[0][0])
-        check = check and (shp_bBox[1][0] > b_bBox[1][0])
-        check = check and (shp_bBox[0][1] > b_bBox[0][1])
-        check = check and (shp_bBox[1][1] < b_bBox[1][1])
+        shp_bBox = [self.tifTrans.get_xy(x_min,y_min),  self.tifTrans.get_xy(x_max,y_max)]
+        check = (shp_bBox[0][0] >= b_bBox[0][0])
+        check = check and (shp_bBox[1][0] <= b_bBox[1][0])
+        check = check and (shp_bBox[0][1] <= b_bBox[0][1])
+        check = check and (shp_bBox[1][1] >= b_bBox[1][1])
 
         if not check:
             raise BoundaryNotInTifError()
@@ -212,8 +234,23 @@ def get_int_box(self, bBox: BBox) -> BBoxInt:
         if not tif_poly.contains(b_poly):
             raise BoundaryNotInTifError()
 
-        return(((x_min,y_min),(x_max,y_max)))
+        return(((x_min, y_min),(x_max, y_max)))
+
 
+    def get_bBox_wgs_84(self, bBox: BBox) -> BBox:
+        """takes a bounding area gets the coordinates of the extremities
+        as if they were clipped by that bounding area
+
+        Args:
+            bBox (BBox): bounding box area to clip within (wgs_84)
+
+        Returns: 
+            BBox: in wgs_84
+        """
+        b = self.get_int_box(bBox)
+        left_top = self.get_wgs_84_coords(b[0][0], b[0][1])
+        right_bottom = self.get_wgs_84_coords(b[1][0], b[1][1])
+        return((left_top, right_bottom))
 
     def read(self) -> np.ndarray:
         """Reade the contents of the geotiff to a zarr array

setup.py

@@ -3,7 +3,7 @@
 import sys
 from setuptools.command.install import install
 
-VERSION = "0.1.3"
+VERSION = "0.1.4"
 
 # Send to pypi
 # python3 setup.py sdist bdist_wheel

tests/test_geotiff.py

@@ -3,36 +3,57 @@
 import os
 from geotiff import GeoTiff
 
+
 @pytest.fixture
 def tiff_file():
     filename = "dem.tif"
     dir = dir = "./tests/inputs/"
     return(os.path.join(dir, filename))
 
+
 @pytest.fixture
-def bounding_box():
+def area_box():
     return([(138.632071411, -32.447310785), (138.644218874, -32.456979174)])
 
 
 @pytest.fixture
 def geoTiff(tiff_file):
     return(GeoTiff(tiff_file))
 
-def test_read(tiff_file, bounding_box, geoTiff: GeoTiff):
+
+def test_read(tiff_file, area_box, geoTiff: GeoTiff):
     print("testing read tiff")
     print(f"reading: {tiff_file}")
-    print(f"Using bBox: {bounding_box}")
-    array = geoTiff.read_box(bounding_box)
+    print(f"Using bBox: {area_box}")
+    array = geoTiff.read_box(area_box)
     print("Sample array:")
     print(array)
     print(array.shape)
     assert isinstance(array, np.ndarray)
 
 
-def test_int_box(bounding_box, geoTiff: GeoTiff):
-    intBox = geoTiff.get_int_box(bounding_box)
+def test_int_box(area_box, geoTiff: GeoTiff):
+    intBox = geoTiff.get_int_box(area_box)
     assert isinstance(intBox, tuple)
     assert len(intBox) == 2
     assert isinstance(intBox[0], tuple)
     assert isinstance(intBox[1], tuple)
-    assert ((125, 143), (169, 178)) == intBox
+    assert ((126, 144), (169, 178)) == intBox
+
+
+def test_conversions(area_box, geoTiff: GeoTiff):
+    int_box = geoTiff.get_int_box(area_box)
+    geoTiff.get_bBox_wgs_84(area_box)
+    i=int_box[0][0] + 5
+    j=int_box[0][1] + 6
+    geoTiff.get_wgs_84_coords(i, j)
+    print(area_box)
+    bounding_box = geoTiff.get_bBox_wgs_84(area_box)
+    assert area_box[0][0] <= bounding_box[0][0]
+    assert area_box[0][1] >= bounding_box[0][1]
+    assert area_box[1][0] >= bounding_box[1][0]
+    assert area_box[1][1] <= bounding_box[1][1]
+
+    print(geoTiff.read_box(area_box).shape)
+    assert geoTiff.tif_bBox_wgs_84 == (
+        (138.5972222222219, -32.40749999999997), (138.69055555555522, -32.49138888888886))