diff --git a/README.md b/README.md index f5073b8..56f7772 100644 --- a/README.md +++ b/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 diff --git a/example.py b/example.py index 9ecaaf0..c6cb0bf 100644 --- a/example.py +++ b/example.py @@ -4,9 +4,11 @@ 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)] @@ -14,9 +16,13 @@ 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") diff --git a/geotiff/__init__.py b/geotiff/__init__.py index 6b5de44..fa2627c 100644 --- a/geotiff/__init__.py +++ b/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 diff --git a/setup.py b/setup.py index 10995e6..699ae5c 100644 --- a/setup.py +++ b/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 diff --git a/tests/test_geotiff.py b/tests/test_geotiff.py index 97ffe91..1f40014 100644 --- a/tests/test_geotiff.py +++ b/tests/test_geotiff.py @@ -3,14 +3,16 @@ 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)]) @@ -18,21 +20,40 @@ def bounding_box(): 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 \ No newline at end of file + 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))