Data extraction#
Extraction pulls values out of a raster, either by narrowing the array to a region of interest (subsetting and clipping) or by sampling pixel values at vector locations (coordinates, points, and polygons). GeoWombat handles the bookkeeping that these tasks usually require, such as translating map coordinates to array indices and reprojecting vector geometries to the raster CRS on-the-fly. The sections below work through each case in turn.
import geowombat as gw
from geowombat.data import rgbn
Subsetting rasters#
Either a rasterio.window.Window
object or tuple can be used with geowombat.open().
Slice a subset using a rasterio.window.Window.
from rasterio.windows import Window
w = Window(row_off=0, col_off=0, height=100, width=100)
bounds = (793475.76, 2049033.03, 794222.03, 2049527.24)
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)
Slice a subset using a tuple of bounded coordinates.
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 Configuration manager 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.
with gw.config.update(ref_bounds=bounds, ref_tar=rgbn):
with gw.open(rgbn) as src:
print(src)
Clipping to bounds#
GeoWombat’s geowombat.clip_by_polygon() is an alternative method to geowombat.config.update. The
geowombat.clip_by_polygon() method limits the bounds of the image to match a polygon, where the polygon
can be a geopandas.GeoDataFrame,
or a path to a file readable with geopandas.read_file.
You can augment the clip by using the argument query on the polygon attributes, and if multiple polygons
are present you can use mask_data to fill nans where polygons are not present, or expand the clip
array bounds by setting expand_by=<n pixels> on each side.
import geowombat as gw
from geowombat.data import l8_224078_20200518, l8_224078_20200518_polygons
import geopandas as gpd
polys = gpd.read_file(l8_224078_20200518_polygons)
with gw.open(l8_224078_20200518) as src:
print(src)
clipped = src.gw.clip_by_polygon(
df,
query="name == water",
mask_data=True,
expand_by=1
)
print(clipped)
Extracting data with coordinates#
To extract values at a coordinate pair, translate the coordinates into array indices.
In [1]: import geowombat as gw
In [2]: from geowombat.data import l8_224078_20200518
# Coordinates in map projection units
In [3]: y, x = -2823031.15, 761592.60
In [4]: with gw.open(l8_224078_20200518) as src:
...: j, i = gw.coords_to_indices(x, y, src)
...: data = src[:, i, j].data.compute()
...:
In [5]: print(data.flatten())
[7448 6882 6090]
A latitude/longitude pair can be extracted after converting to the map projection.
In [6]: import geowombat as gw
In [7]: from geowombat.data import l8_224078_20200518
# Coordinates in latitude/longitude
In [8]: lat, lon = -25.50142964, -54.39756038
In [9]: with gw.open(l8_224078_20200518) as src:
...: x, y = gw.lonlat_to_xy(lon, lat, src)
...: j, i = gw.coords_to_indices(x, y, src)
...: data = src[:, i, j].data.compute()
...:
In [10]: 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().
In [11]: import geowombat as gw
In [12]: from geowombat.data import l8_224078_20200518, l8_224078_20200518_points
In [13]: with gw.open(l8_224078_20200518) as src:
....: df = src.gw.extract(l8_224078_20200518_points)
....:
In [14]: 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.23) 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 would have transformed the CRS on-the-fly.
In [15]: import geowombat as gw
In [16]: from geowombat.data import l8_224078_20200518, l8_224078_20200518_points
In [17]: import geopandas as gpd
In [18]: point_df = gpd.read_file(l8_224078_20200518_points)
In [19]: print(point_df.crs)
EPSG:32621
# Transform the CRS to WGS84 lat/lon
In [20]: point_df = point_df.to_crs('epsg:4326')
In [21]: print(point_df.crs)
epsg:4326
In [22]: with gw.open(l8_224078_20200518) as src:
....: df = src.gw.extract(point_df)
....:
In [23]: 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.23) 5 7398 6711 6007
Set the data band names.
In [24]: import geowombat as gw
In [25]: from geowombat.data import l8_224078_20200518, l8_224078_20200518_points
In [26]: 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()
....: )
....:
In [27]: 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.23) 5 7398 6711 6007
Extracting data with polygon geometry#
To extract values within polygons, use the same geowombat.extract() function.
In [28]: from geowombat.data import l8_224078_20200518, l8_224078_20200518_polygons
In [29]: 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()
....: )
....:
In [30]: print(df)
id point geometry name blue green red
0 0 0 POINT (737559.502 -2795247.772) water 7994 7423 6272
1 0 1 POINT (737589.502 -2795247.772) water 8017 7428 6292
2 0 2 POINT (737619.502 -2795247.772) water 8008 7446 6292
3 0 3 POINT (737649.502 -2795247.772) water 8008 7412 6291
4 0 4 POINT (737679.502 -2795247.772) water 8018 7398 6250
.. .. ... ... ... ... ... ...
667 3 667 POINT (739038.667 -2811819.677) developed 8567 8564 8447
668 3 668 POINT (739068.667 -2811819.677) developed 8099 7676 7332
669 3 669 POINT (739098.667 -2811819.677) developed 10151 9651 10153
670 3 670 POINT (739128.667 -2811819.677) developed 8065 7735 7501
671 3 671 POINT (739158.667 -2811819.677) developed 9343 8987 9247
[672 rows x 7 columns]