Skip to content

GeoTensor numpy API

  • Author: Gonzalo Mateo-García

This tutorial shows how to operate with the numpy API with GeoTensor objects. The GeoTensor class is a subclass of np.ndarray that stores the spatial affine transform, the coordinate reference system and no data value (called fill_value_default). When a GeoTensor is sliced, its transform attribute is shifted accordingly. Additionally its transform is also shifted if the GeoTensor is padded. GeoTensors are restricted to be 2D, 3D or 4D arrays and their two last dimensions are assumed to be the y and x spatial axis.

As a subclass of np.ndarray, operations with GeoTensor objects work similar than operations with np.ndarrays. However, there are some restrictions that we have implemented to keep consistency with the GeoTensor concept. If you need to use the numpy implementation you can access the bare numpy object with the .values attribute. Below there's a list with restrictions on numpy operations:

  1. Slicing a GeoTensor is more restrictive than a numpy array. It only allows to slice with lists, numbers or slices. In particular the spatial dimensions can only be sliced with slices. Slicing for inplace modification is not restricted (i.e. you can slice with boolean arrays to modify certain values of the object).
  2. Binary operations (such as add, multiply etc) check for GeoTensor inputs if they have the same_extent; that is, same transform crs and spatial dimensions (width and height).
  3. squeeze, expand_dims and transpose make sure spatial dimensions (last two axes) are not modified and kept at the end of the array.
  4. concatenate and stack make sure all operated GeoTensors have same_extent and shape. concatenate does not allow to concatenate on the spatial dims.
  5. Reductions (such as np.mean or np.all) return GeoTensor object if the spatial dimensions are preserved and np.ndarray or scalars otherwise.

Load some sample data

from georeader import plot
from georeader.geotensor import GeoTensor
from georeader import save
import os
import numpy as np

s2url = "https://huggingface.co/datasets/isp-uv-es/WorldFloodsv2/resolve/main/test/S2/EMSR264_18MIANDRIVAZODETAIL_DEL_v2.tif"
filelocal = os.path.basename(s2url)
if not os.path.exists(filelocal):
    s2img = GeoTensor.load_file(s2url)
    save.save_tiled_geotiff(s2img, filelocal)
else:
    s2img = GeoTensor.load_file(filelocal)

# Set nodata value to zero
s2img.fill_value_default = 0
s2img

 
         Transform: | 10.00, 0.00, 537430.00|
| 0.00,-10.00, 7844180.00|
| 0.00, 0.00, 1.00|
         Shape: (15, 1313, 1530)
         Resolution: (10.0, 10.0)
         Bounds: (537430.0, 7831050.0, 552730.0, 7844180.0)
         CRS: EPSG:32738
         fill_value_default: 0

Padding and slicing as with numpy arrays

transform_before = s2img.transform
s2img = s2img.pad(((0,0), (25, 25), (10, 10)))
s2img

 
         Transform: | 10.00, 0.00, 537330.00|
| 0.00,-10.00, 7844430.00|
| 0.00, 0.00, 1.00|
         Shape: (15, 1363, 1550)
         Resolution: (10.0, 10.0)
         Bounds: (537330.0, 7830800.0, 552830.0, 7844430.0)
         CRS: EPSG:32738
         fill_value_default: 0

Padding changes transform

assert s2img.transform != transform_before, "adding changes transform"
s2img.transform, transform_before

(Affine(10.0, 0.0, 537330.0,
        0.0, -10.0, 7844430.0),
 Affine(10.0, 0.0, 537430.0,
        0.0, -10.0, 7844180.0))

rgb = (s2img[[3,2,1]] / 3_500).clip(0,1)
plot.show(rgb)

<Axes: >
No description has been provided for this image

validmask and invalidmask methods

invalids = np.all(s2img.invalidmask(),axis=0) 

assert isinstance(invalids, GeoTensor), "invalids is of class GeoTensor!"

plot.plot_segmentation_mask(invalids, interpretation_array=["valid", "invalid"],
                           color_array=["C0","C1"])

<Axes: >
No description has been provided for this image

Slicing propagates the affine transform

rgb_slice = rgb[:, 20:300, 50:340]
rgb_slice

 
         Transform: | 10.00, 0.00, 537830.00|
| 0.00,-10.00, 7844230.00|
| 0.00, 0.00, 1.00|
         Shape: (3, 280, 290)
         Resolution: (10.0, 10.0)
         Bounds: (537830.0, 7841430.0, 540730.0, 7844230.0)
         CRS: EPSG:32738
         fill_value_default: 0
        

assert rgb_slice.transform != rgb.transform, f"Different transform"

rgb_slice.transform, rgb.transform

(Affine(10.0, 0.0, 537830.0,
        0.0, -10.0, 7844230.0),
 Affine(10.0, 0.0, 537330.0,
        0.0, -10.0, 7844430.0))

Slicing with step

rgb_slice_2 = rgb[:, slice(20,300,2), slice(50,340)]
rgb_slice_2

 
         Transform: | 10.00, 0.00, 537830.00|
| 0.00,-20.00, 7844230.00|
| 0.00, 0.00, 1.00|
         Shape: (3, 140, 290)
         Resolution: (10.0, 20.0)
         Bounds: (537830.0, 7841430.0, 540730.0, 7844230.0)
         CRS: EPSG:32738
         fill_value_default: 0
        

plot.show(rgb_slice_2)

<Axes: >
No description has been provided for this image

Reversing

rgb_slice_rev = rgb[:, slice(300,20,-1), slice(50,340)]
rgb_slice_rev

 
         Transform: | 10.00, 0.00, 537830.00|
| 0.00, 10.00, 7841430.00|
| 0.00, 0.00, 1.00|
         Shape: (3, 280, 290)
         Resolution: (10.0, 10.0)
         Bounds: (537830.0, 7841430.0, 540730.0, 7844230.0)
         CRS: EPSG:32738
         fill_value_default: 0
        

plot.show(rgb_slice_rev)

<Axes: >
No description has been provided for this image

np.mean and other reduction ops return GeoTensor

When reduction axis are not the spatial axis, the returned object is a GeoTensor

import numpy as np

greyimg = np.mean(rgb, axis=0)
greyimg

 
         Transform: | 10.00, 0.00, 537330.00|
| 0.00,-10.00, 7844430.00|
| 0.00, 0.00, 1.00|
         Shape: (1363, 1550)
         Resolution: (10.0, 10.0)
         Bounds: (537330.0, 7830800.0, 552830.0, 7844430.0)
         CRS: EPSG:32738
         fill_value_default: 0
        

plot.show(greyimg, add_colorbar_next_to=True)

<Axes: >
No description has been provided for this image

When reducing over the spatial axis, return object is np.ndarray:

average_reflectance_band = np.mean(rgb, axis=(-2,-1))
average_reflectance_band

array([0.39475575, 0.36224817, 0.34591881])

vectorized ufuncs

greyimg_sin = np.sin(greyimg)
greyimg_sin

 
         Transform: | 10.00, 0.00, 537330.00|
| 0.00,-10.00, 7844430.00|
| 0.00, 0.00, 1.00|
         Shape: (1363, 1550)
         Resolution: (10.0, 10.0)
         Bounds: (537330.0, 7830800.0, 552830.0, 7844430.0)
         CRS: EPSG:32738
         fill_value_default: 0

Ops of GeoTensor with numbers and np.ndarray

clouds = greyimg &gt;= .9
plot.plot_segmentation_mask(clouds, interpretation_array=["clear","cloud"], 
                            color_array=["C0", "white"])

<Axes: >
No description has been provided for this image 
cloud_percentage = np.mean(clouds)
cloud_percentage

np.float64(0.022976830047570587)

Example: normalize by bands

import matplotlib.pyplot as plt

average_reflectance_band = np.mean(rgb, axis=(-2,-1), keepdims=True)
std_reflectance = np.mean(rgb, axis=(-2,-1), keepdims=True)

rgb_norm = (rgb - average_reflectance_band) / std_reflectance

fig, ax = plt.subplots(1, 3, figsize=(15, 5),sharey=True, tight_layout=True)
for i in range(3):
    plot.show(rgb_norm[i], add_colorbar_next_to=True,ax=ax[i], title=f"B{i+2}")
No description has been provided for this image

Computing RS indexes

b11 = s2img[11].astype(np.float64)
b3 = s2img[2].astype(np.float64)

mndwi = (b11 - b3) / (b11 + b3)
plot.show(mndwi, add_colorbar_next_to=True)

/home/gonzalo/git/georeader/georeader/geotensor.py:556: RuntimeWarning: invalid value encountered in divide
  result_values = self.values / other


<Axes: >
No description has been provided for this image 
water = mndwi &lt; 0
land_water_clouds = GeoTensor(np.ones(clouds.shape, dtype=np.uint8),
                              fill_value_default=0,
                              crs=clouds.crs,
                              transform=clouds.transform)

land_water_clouds[water] = 2
land_water_clouds[clouds] = 3
land_water_clouds[invalids] = 0

plot.plot_segmentation_mask(land_water_clouds, 
                            interpretation_array=["invalids","clear","water","cloud"], 
                            color_array=["#000000","#c66d43","#437cc6","#eeeeee"])

<Axes: >
No description has been provided for this image

.values attribute returns a view of the object as np.ndarray

clouds.values

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])

Apply function that expects np.ndarray objects

We can pass a GeoTensor to any function that expects np.ndarray. Depending on the function the return type may be a GeoTensor or a np.ndarray. If the return type is np.ndarray we can wrap the array with array_as_geotensor method to convert back to GeoTensor object with the same spatial info.

from scipy import ndimage as ndi

greyimg_smoothed = ndi.gaussian_filter(greyimg, sigma=5, cval=0, mode="reflect")
greyimg_smoothed

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

greyimg_smoothed = rgb.array_as_geotensor(greyimg_smoothed)
plot.show(greyimg_smoothed, add_colorbar_next_to=True)

<Axes: >
No description has been provided for this image

Licence

The georeader package is published under a GNU Lesser GPL v3 licence

georeader tutorials and notebooks are released under a Creative Commons non-commercial licence.

If you find this work useful please cite:

@article{portales-julia_global_2023,
    title = {Global flood extent segmentation in optical satellite images},
    volume = {13},
    issn = {2045-2322},
    doi = {10.1038/s41598-023-47595-7},
    number = {1},
    urldate = {2023-11-30},
    journal = {Scientific Reports},
    author = {Portalés-Julià, Enrique and Mateo-García, Gonzalo and Purcell, Cormac and Gómez-Chova, Luis},
    month = nov,
    year = {2023},
    pages = {20316},
}