Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,11 @@ test-cov: ## Run tests with coverage report
@echo "🧪 Running tests with coverage"
@poetry run pytest tests/ -v --cov=georeader --cov-report=term-missing --cov-report=html

.PHONY: test-doctest
test-doctest: ## Run doctests in core modules
@echo "🧪 Running doctests"
@poetry run pytest --doctest-modules georeader/geotensor.py -v

.PHONY: build
build: clean-build ## Build wheel file using poetry
@echo "📦 Building wheel file"
Expand Down
98 changes: 62 additions & 36 deletions georeader/geotensor.py
Original file line number Diff line number Diff line change
Expand Up @@ -651,7 +651,7 @@ def __add__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> # Create sample GeoTensor (3 bands, 100x100 pixels)
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
Expand All @@ -660,7 +660,8 @@ def __add__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
>>>
>>> # Add scalar to all pixels
>>> gt_offset = gt + 0.1
>>> print(gt_offset.shape) # (3, 100, 100)
>>> print(gt_offset.shape)
(3, 100, 100)
>>>
>>> # Add per-band offset using broadcasting
>>> band_offsets = np.array([0.1, 0.2, 0.3])[:, None, None] # Shape: (3, 1, 1)
Expand Down Expand Up @@ -726,7 +727,7 @@ def __sub__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> # Create sample GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
Expand Down Expand Up @@ -791,7 +792,7 @@ def __mul__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100)
Expand Down Expand Up @@ -866,7 +867,7 @@ def __truediv__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100)
Expand Down Expand Up @@ -1241,6 +1242,11 @@ def expand_dims(self, axis: Union[int, tuple]) -> Self:
ValueError: If trying to add dimensions at or after the spatial dimensions.

Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader.geotensor import GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> crs = "EPSG:32630"
>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> # Add a new dimension at axis 0
>>> gt_expanded = gt.expand_dims(0)
Expand Down Expand Up @@ -1375,14 +1381,14 @@ def isel(self, sel: Dict[str, Union[slice, list, int]]) -> Self:
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> # Create 3-band 100x100 GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100).astype(np.float32)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>> print(f"Original: {gt.shape}, origin: ({gt.transform.c}, {gt.transform.f})")
Original: (3, 100, 100), origin: (500000, 4650000)
Original: (3, 100, 100), origin: (500000.0, 4650000.0)
>>>
>>> # Slice spatial subset (maintains bands)
>>> subset = gt.isel({"x": slice(20, 80), "y": slice(10, 60)})
Expand Down Expand Up @@ -1492,33 +1498,29 @@ def footprint(self, crs: Optional[str] = None) -> Polygon:

Examples:
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>> import numpy as np
>>>
>>> # Example 1: Get footprint in native CRS
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000) # UTM transform
>>> data = np.random.rand(3, 100, 100)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>> footprint_utm = gt.footprint()
>>> print(f"Bounds (UTM): {footprint_utm.bounds}")
>>> # (xmin, ymin, xmax, ymax) in UTM meters: (500000, 4649000, 501000, 4650000)
>>> footprint_utm.bounds
(500000.0, 4649000.0, 501000.0, 4650000.0)

>>> # Example 2: Transform footprint to WGS84 for web mapping
>>> footprint_wgs84 = gt.footprint(crs="EPSG:4326")
>>> print(f"Bounds (WGS84): {footprint_wgs84.bounds}")
>>> # (lon_min, lat_min, lon_max, lat_max) in degrees
>>> # Bounds returned in degrees: (lon_min, lat_min, lon_max, lat_max)
>>> # Can be used directly in Leaflet, Google Maps, etc.

>>> # Example 3: Check if rasters overlap
>>> gt1 = GeoTensor.load_file('image1.tif')
>>> gt2 = GeoTensor.load_file('image2.tif')
>>> # Get footprints in common CRS
>>> fp1 = gt1.footprint(crs="EPSG:4326")
>>> fp2 = gt2.footprint(crs="EPSG:4326")
>>> if fp1.intersects(fp2):
... print("Rasters overlap!")
>>> gt1 = GeoTensor.load_file('image1.tif') # doctest: +SKIP
>>> gt2 = GeoTensor.load_file('image2.tif') # doctest: +SKIP
>>> fp1 = gt1.footprint(crs="EPSG:4326") # doctest: +SKIP
>>> fp2 = gt2.footprint(crs="EPSG:4326") # doctest: +SKIP
>>> if fp1.intersects(fp2): # doctest: +SKIP
... overlap_area = fp1.intersection(fp2).area
... print(f"Overlap area: {overlap_area} square degrees")

>>> # Example 4: Export footprint as GeoJSON
>>> from shapely.geometry import mapping
Expand All @@ -1529,21 +1531,20 @@ def footprint(self, crs: Optional[str] = None) -> Polygon:
... "geometry": mapping(footprint),
... "properties": {"name": "Raster extent"}
... }
>>> with open('footprint.geojson', 'w') as f:
>>> with open('footprint.geojson', 'w') as f: # doctest: +SKIP
... json.dump(geojson, f)

>>> # Example 5: Check if point is within raster extent
>>> from shapely.geometry import Point
>>> point_of_interest = Point(-3.7038, 40.4168) # Madrid coordinates
>>> footprint = gt.footprint(crs="EPSG:4326")
>>> if footprint.contains(point_of_interest):
... print("Point is within raster extent")
>>> contains_point = footprint.contains(point_of_interest)

>>> # Example 6: Calculate raster area in square kilometers
>>> footprint_utm = gt.footprint() # In UTM (meters)
>>> area_sqm = footprint_utm.area
>>> area_sqkm = area_sqm / 1_000_000
>>> print(f"Raster covers {area_sqkm:.2f} km²")
>>> area_sqkm = footprint_utm.area / 1_000_000
>>> round(area_sqkm, 2)
1.0

Note:
- Footprint represents the full raster extent, including nodata regions
Expand Down Expand Up @@ -1696,7 +1697,7 @@ def pad(
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> # Create 3-band 100x100 GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
Expand Down Expand Up @@ -1767,9 +1768,14 @@ def pad_array(
GeoTensor: padded GeoTensor.

Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader.geotensor import GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> crs = "EPSG:32630"
>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt.pad_array({"x": (10, 10), "y": (10, 10)})
>>> assert gt.shape == (3, 120, 120)
>>> padded = gt.pad_array({"x": (10, 10), "y": (10, 10)})
>>> assert padded.shape == (3, 120, 120)
"""
if constant_values is None and mode == "constant":
if self.fill_value_default is None:
Expand Down Expand Up @@ -1842,10 +1848,15 @@ def resize(
resized GeoTensor

Examples:
>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> resized = gt.resize((50, 50))
>>> assert resized.shape == (3, 50, 50)
>>> assert resized.res == (2*gt.res[0], 2*gt.res[1])
>>> import numpy as np # doctest: +SKIP
>>> import rasterio # doctest: +SKIP
>>> from georeader.geotensor import GeoTensor # doctest: +SKIP
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000) # doctest: +SKIP
>>> crs = "EPSG:32630" # doctest: +SKIP
>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs) # doctest: +SKIP
>>> resized = gt.resize((50, 50)) # doctest: +SKIP
>>> assert resized.shape == (3, 50, 50) # doctest: +SKIP
>>> assert resized.res == (2*gt.res[0], 2*gt.res[1]) # doctest: +SKIP
"""
input_shape = self.shape
spatial_shape = input_shape[-2:]
Expand Down Expand Up @@ -1956,6 +1967,11 @@ def transpose(self, axes=None) -> Self:
ValueError: If the spatial dimensions are moved from their last positions.

Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader.geotensor import GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> crs = "EPSG:32630"
>>> gt = GeoTensor(np.random.rand(3, 4, 100, 200), transform, crs)
>>> # Shape is (3, 4, 100, 200)
>>> gt_t = gt.transpose()
Expand Down Expand Up @@ -2189,6 +2205,11 @@ def write_from_window(self, data: np.ndarray, window: rasterio.windows.Window):
window: Window object that specifies the spatial location to write the data

Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader.geotensor import GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> crs = "EPSG:32630"
>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> data = np.random.rand(3, 50, 50)
>>> window = rasterio.windows.Window(col_off=7, row_off=9, width=50, height=50)
Expand Down Expand Up @@ -2273,7 +2294,7 @@ def read_from_window(
>>> import numpy as np
>>> import rasterio
>>> import rasterio.windows
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> # Create test GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
Expand Down Expand Up @@ -2403,7 +2424,7 @@ def stack(cls, geotensors: List[Self]) -> Self:
Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>> from georeader.geotensor import GeoTensor
>>>
>>> # Create multiple GeoTensors (e.g., different dates)
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
Expand Down Expand Up @@ -2490,10 +2511,15 @@ def concatenate(cls, geotensors: List[Self], axis: int = 0) -> Self:
geotensor with extra dim at the front: (len(geotensors),) + shape

Examples:
>>> import numpy as np
>>> import rasterio
>>> from georeader.geotensor import GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> crs = "EPSG:32630"
>>> gt1 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt3 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt = concatenate([gt1, gt2, gt3], axis=0)
>>> gt = GeoTensor.concatenate([gt1, gt2, gt3], axis=0)
>>> assert gt.shape == (9, 100, 100)
"""
assert len(geotensors) > 0, "Empty list provided can't concat"
Expand Down
1 change: 1 addition & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@ show_error_codes = "True"

[tool.pytest.ini_options]
testpaths = ["tests"]
doctest_optionflags = ["ELLIPSIS", "NORMALIZE_WHITESPACE"]

[tool.ruff]
target-version = "py39"
Expand Down