Learning Objectives

  • Subset bands by index or name

  • Extract raster data by row and column number

  • Extract data by bounding window

  • Extract raster data by coordinates

  • Extract raster data by geometry (point, polygon)


Raster Data Extraction#

Raster data is often of little use unless we can extract and summarize the data. For instance, extracting raster to points by coordinates allows us to pass data to machine learning models for land cover classification or cloud masking.

Subsetting rasters#

We can subset sections of the data array in a number of ways. In this case we will create a slice based on row and column location to subset LandSat data using a rasterio.window.Window.

Either a rasterio.window.Window object or tuple can be used with geowombat.open.

import geowombat as gw
from geowombat.data import rgbn

from rasterio.windows import Window
w = Window(row_off=0, col_off=0, height=100, width=100)

with gw.open(rgbn,
                band_names=['blue', 'green', 'red'],
                num_workers=8,
                indexes=[1, 2, 3],
                window=w,
                out_dtype='float32') as src:
    print(src)
<xarray.DataArray 'array-c1686f70048cb766e2fe94fb3c34fb6d' (band: 3, y: 100, x: 100)>
dask.array<array, shape=(3, 100, 100), dtype=float32, chunksize=(1, 64, 64), chunktype=numpy.ndarray>
Coordinates:
  * band     (band) <U5 'blue' 'green' 'red'
  * y        (y) float64 2.05e+06 2.05e+06 2.05e+06 ... 2.05e+06 2.05e+06
  * x        (x) float64 7.93e+05 7.93e+05 7.93e+05 ... 7.935e+05 7.935e+05
Attributes:
    transform:          | 5.00, 0.00, 792988.00|\n| 0.00,-5.00, 2050382.00|\n...
    crs:                +init=epsg:32618
    res:                (5.0, 5.0)
    is_tiled:           1
    nodatavals:         (nan, nan, nan, nan)
    offsets:            (0.0, 0.0, 0.0, 0.0)
    data_are_separate:  0
    data_are_stacked:   0

We can also slice a subset of data using a tuple of bounded coordinates.

bounds = (793475.76, 2049033.03, 794222.03, 2049527.24)

with gw.open(rgbn,
                band_names=['green', 'red', 'nir'],
                num_workers=8,
                indexes=[2, 3, 4],
                bounds=bounds,
                out_dtype='float32') as src:
    print(src)

The configuration manager provides an alternative method to subset rasters. See tutorial-config for more details.

with gw.config.update(ref_bounds=bounds):

    with gw.open(rgbn) as src:
        print(src)

By default, the subset will be returned by the upper left coordinates of the bounds, potentially shifting cell alignment with the reference raster. To subset a raster and align it to the same grid, use the ref_tar keyword. This is equivalent to a “snap raster” in ArcGIS.

with gw.config.update(ref_bounds=bounds, ref_tar=rgbn):

    with gw.open(rgbn) as src:
        print(src)

Extracting data by coordinates#

To extract values at a coordinate pair, translate the coordinates into array indices. For extraction by geometry, for instance with a shapefile, see extract by point geometry.

import geowombat as gw
from geowombat.data import l8_224078_20200518

# Coordinates in map projection units
y, x = -2823031.15, 761592.60

with gw.open(l8_224078_20200518) as src:
    # Transform the map coordinates to data indices
    j, i = gw.coords_to_indices(x, y, src)
    # Subset by index
    data = src[:, i, j].data.compute()

print(data.flatten())
[7448 6882 6090]

A latitude/longitude pair can be extracted after converting to the map projection.

import geowombat as gw
from geowombat.data import l8_224078_20200518

# Coordinates in latitude/longitude
lat, lon = -25.50142964, -54.39756038

with gw.open(l8_224078_20200518) as src:
    # Transform the coordinates to map units
    x, y = gw.lonlat_to_xy(lon, lat, src)
    # Transform the map coordinates to data indices
    j, i = gw.coords_to_indices(x, y, src)
    data = src[:, i, j].data.compute()

print(data.flatten())
[7448 6882 6090]

Extracting data with point geometry#

In the example below, ‘l8_224078_20200518_points’ is a GeoPackage of point locations, and the output df is a GeoPandas GeoDataFrame. To extract the raster values at the point locations, use geowombat.extract.

import geowombat as gw
from geowombat.data import l8_224078_20200518, l8_224078_20200518_points

with gw.open(l8_224078_20200518) as src:
    df = src.gw.extract(l8_224078_20200518_points)

print(df)
        name                         geometry  id     1     2     3
0      water  POINT (741522.314 -2811204.698)   0  7966  7326  6254
1       crop  POINT (736140.845 -2806478.364)   1  8030  7490  8080
2       tree  POINT (745919.508 -2805168.579)   2  7561  6874  6106
3  developed  POINT (739056.735 -2811710.662)   3  8302  8202  8111
4      water  POINT (737802.183 -2818016.412)   4  8277  7982  7341
5       tree  POINT (759209.443 -2828566.230)   5  7398  6711  6007

Note

The line df = src.gw.extract(l8_224078_20200518_points) could also have been written as df = gw.extract(src, l8_224078_20200518_points).

In the previous example, the point vector had a CRS that matched the raster (i.e., EPSG=32621, or UTM zone 21N). If the CRS had not matched, the geowombat.extract function transforms the CRS on-the-fly.

import geowombat as gw
from geowombat.data import l8_224078_20200518, l8_224078_20200518_points
import geopandas as gpd

point_df = gpd.read_file(l8_224078_20200518_points)
print(point_df.crs)

# Transform the CRS to WGS84 lat/lon
point_df = point_df.to_crs('epsg:4326')
print(point_df.crs)

with gw.open(l8_224078_20200518) as src:
    df = src.gw.extract(point_df)

print(df)
epsg:32621
epsg:4326
        name                         geometry  id     1     2     3
0      water  POINT (741522.314 -2811204.698)   0  7966  7326  6254
1       crop  POINT (736140.845 -2806478.364)   1  8030  7490  8080
2       tree  POINT (745919.508 -2805168.579)   2  7561  6874  6106
3  developed  POINT (739056.735 -2811710.662)   3  8302  8202  8111
4      water  POINT (737802.183 -2818016.412)   4  8277  7982  7341
5       tree  POINT (759209.443 -2828566.230)   5  7398  6711  6007

Set the data band names using sensor = 'bgr', which assigns the band names blue, green, red.

import geowombat as gw
from geowombat.data import l8_224078_20200518, l8_224078_20200518_points

with gw.config.update(sensor='bgr'):
    with gw.open(l8_224078_20200518) as src:
        df = src.gw.extract(l8_224078_20200518_points,
                            band_names=src.band.values.tolist())

print(df)
        name                         geometry  id  blue  green   red
0      water  POINT (741522.314 -2811204.698)   0  7966   7326  6254
1       crop  POINT (736140.845 -2806478.364)   1  8030   7490  8080
2       tree  POINT (745919.508 -2805168.579)   2  7561   6874  6106
3  developed  POINT (739056.735 -2811710.662)   3  8302   8202  8111
4      water  POINT (737802.183 -2818016.412)   4  8277   7982  7341
5       tree  POINT (759209.443 -2828566.230)   5  7398   6711  6007

Extracting time series images by point geometry#

We can also easily extract a time series of raster images. Extracted pixel values are provided in ‘wide’ format with appropriate labels, for instance the column ‘t2_blue’ would be the blue band for the second time period

from geowombat.data import l8_224078_20200518, l8_224078_20200518_points

with gw.config.update(sensor='bgr'):
    with gw.open([l8_224078_20200518, l8_224078_20200518],
            time_names=['t1', 't2'],
            stack_dim='time') as src:

        # Extract and by point geometry
        df = src.gw.extract(l8_224078_20200518_points)

print(df)
        name                         geometry  id  t1_blue  t1_green  t1_red  \
0      water  POINT (741522.314 -2811204.698)   0     7966      7326    6254   
1       crop  POINT (736140.845 -2806478.364)   1     8030      7490    8080   
2       tree  POINT (745919.508 -2805168.579)   2     7561      6874    6106   
3  developed  POINT (739056.735 -2811710.662)   3     8302      8202    8111   
4      water  POINT (737802.183 -2818016.412)   4     8277      7982    7341   
5       tree  POINT (759209.443 -2828566.230)   5     7398      6711    6007   

   t2_blue  t2_green  t2_red  
0     7966      7326    6254  
1     8030      7490    8080  
2     7561      6874    6106  
3     8302      8202    8111  
4     8277      7982    7341  
5     7398      6711    6007  

Extracting data by polygon geometry#

To extract values within polygons, use the same geowombat.extract function.

from geowombat.data import l8_224078_20200518, l8_224078_20200518_polygons

with gw.config.update(sensor='bgr'):
    with gw.open(l8_224078_20200518) as src:
        df = src.gw.extract(l8_224078_20200518_polygons,
                            band_names=src.band.values.tolist())

    print(df)
     id  point                         geometry       name   blue  green  \
0     0      0  POINT (737535.000 -2795205.000)      water   8017   7435   
1     0      1  POINT (737565.000 -2795205.000)      water   8016   7439   
2     0      2  POINT (737595.000 -2795205.000)      water   8012   7442   
3     0      3  POINT (737625.000 -2795205.000)      water   7997   7422   
4     0      4  POINT (737655.000 -2795205.000)      water   7997   7405   
..   ..    ...                              ...        ...    ...    ...   
667   3    667  POINT (739005.000 -2811795.000)  developed   9014   8236   
668   3    668  POINT (739035.000 -2811795.000)  developed   8567   8564   
669   3    669  POINT (739065.000 -2811795.000)  developed   8099   7676   
670   3    670  POINT (739095.000 -2811795.000)  developed  10151   9651   
671   3    671  POINT (739125.000 -2811795.000)  developed   8065   7735   

       red  
0     6283  
1     6294  
2     6295  
3     6284  
4     6266  
..     ...  
667   8325  
668   8447  
669   7332  
670  10153  
671   7501  

[672 rows x 7 columns]

Calculate mean pixel value by polygon#

It is simple then to calculate the mean value of pixels within each polygon by using the polygon id column and pandas groupby function. You can easily calculate other statistics like min, max, median etc.

from geowombat.data import l8_224078_20200518, l8_224078_20200518_polygons

with gw.config.update(sensor='bgr'):
    with gw.open(l8_224078_20200518) as src:
        df = src.gw.extract(l8_224078_20200518_polygons,
                            band_names=src.band.values.tolist())
        # use pandas groupby to calc pixel mean  
        df = df.groupby('id').mean()
    print(df)
    point         blue        green          red
id                                              
0   103.5  7990.052885  7388.432692  6264.807692
1   304.0  7692.481865  7037.419689  7571.207254
2   497.0  7506.901554  6838.704663  6091.932642
3   632.5  8668.423077  8294.717949  8312.192308