Learning Objectives

  • Reproject remotely sensed data (change CRS)

  • Reproject on-the-fly

  • Compare resampling options


Configuration manager#

What is a context manager?#

In short, a context manager ensures proper file closing using with statements. But it also allows us to set up default behavoirs for opening and writing our images.

What is the purpose of GeoWombat’s context manager?#

The examples shown in reading remotely sensed data opened the entire raster as DataArrays as they were stored on file. The configuration manager allows easy control over opened raster dimensions, alignment, and transformations.

For instance you might want to set the bound (extent) of your analysis. By setting bounds with the configuration manager we will minimize our overhead (less data processed) and uniformly treat all images we process.

How do I use it?#

To use GeoWombat’s configuration manager, just call geowombat.config.update before opening a file. For example,

import geowombat as gw

with gw.config.update(<keywords>...):

# Every file opened within the configuration block will use
# configuration keywords
with gw.open('image.tif') as src:
    # do something

geowombat.config.update stores keywords in a dictionary. To see all GeoWombat configuration keywords, just iterate over the dictionary.

import geowombat as gw
from geowombat.data import l8_224078_20200518

# Using the manager without keywords will set defaults
with gw.config.update():
    with gw.open(l8_224078_20200518) as src:
        for k, v in src.gw.config.items():
            print('Keyword:', k.ljust(15), 'Value:', v)
Keyword: with_config     Value: True
Keyword: ignore_warnings Value: False
Keyword: sensor          Value: None
Keyword: scale_factor    Value: None
Keyword: offset          Value: None
Keyword: nodata          Value: None
Keyword: ref_image       Value: None
Keyword: ref_bounds      Value: None
Keyword: ref_crs         Value: None
Keyword: ref_res         Value: None
Keyword: ref_tar         Value: None
Keyword: blockxsize      Value: 512
Keyword: blockysize      Value: 512
Keyword: compress        Value: None
Keyword: driver          Value: GTiff
Keyword: tiled           Value: True
Keyword: bigtiff         Value: NO
Keyword: l57_angles_path Value: None
Keyword: l8_angles_path  Value: None
Keyword: nasa_earthdata_user Value: None
Keyword: nasa_earthdata_key Value: None
Keyword: nasa_earthdata_code Value: None

Image Properties#

Certain raster properties can be set via the config.update. This includes setting no data values - to be masked, and pixel value scaling via scale_factor.

# treat 0 as missing value
with gw.config.update(nodata=0):
  with gw.open(l8_224078_20200518) as src:
    print(src.attrs['nodatavals'])
    #use src=src.gw.mask_nodata() to replace 0 with nan
(0, 0, 0)
# multiply pixel values by 0.0001
with gw.config.update(scale_factor=0.0001):
  with gw.open(l8_224078_20200518) as src:
    print(src.attrs['scales'])
(1.0, 1.0, 1.0)

Reference settings: CRS#

Configuration keywords beginning with ref are the most important commands when opening rasters. For example, to transform the CRS of the data on-the-fly, use ref_crs. For more on Coordinate Reference Systems, see here.

import geowombat as gw
from geowombat.data import l8_224078_20200518

proj4 = "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs "

# Without the manager
with gw.open(l8_224078_20200518) as src:
    print(src.crs)

# With the manager
with gw.config.update(ref_crs=proj4):
    with gw.open(l8_224078_20200518) as src:
        print(src.crs)
32621
+proj=aea +lat_0=-32 +lon_0=-60 +lat_1=-5 +lat_2=-42 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs +type=crs

Reference settings: Cell size#

It is possible to combine multiple configuration keywords. In the example below, the raster CRS is transformed from UTM to Albers Equal Area with a resampled cell size of 100m x 100m.

import geowombat as gw
from geowombat.data import l8_224078_20200518

# Without the manager
with gw.open(l8_224078_20200518) as src:
    print(src.gw.celly, src.gw.cellx)

# With the manager
with gw.config.update(ref_crs=proj4, ref_res=(100, 100)):
    with gw.open(l8_224078_20200518) as src:
        print(src.gw.celly, src.gw.cellx)
30.0 30.0
100.0 100.0

Reference settings: Bounds#

To subset an image, specify bounds as a tuple of (left, bottom, right, top) or a rasterio BoundingBox object.

import geowombat as gw
from geowombat.data import l8_224078_20200518
from rasterio.coords import BoundingBox

bounds = BoundingBox(left=724634.17, bottom=-2806501.39, right=737655.48, top=-2796221.42)

# or
# bounds = (724634.17, -2806501.39, 737655.48, -2796221.42)

# Without the manager
with gw.open(l8_224078_20200518) as src:
    print(src.gw.bounds)

# With the manager
with gw.config.update(ref_bounds=bounds):
    with gw.open(l8_224078_20200518) as src:
        print(src.gw.bounds)
(717345.0, -2832795.0, 778575.0, -2776995.0)
(724634.17, -2806481.42, 737654.17, -2796221.42)

Reference settings: Snap Raster Target#

By default, the bounding 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)

Reference Image#

To use another image as a reference, just set ref_image. Then, the opened file’s bounds, CRS, and cell size will be transformed to match those of the reference image.

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

# Without the manager
with gw.open(l8_224078_20200518) as src:
    print(src.gw.bounds)

with gw.open(l8_224077_20200518_B2) as src:
    print(src.gw.bounds)

# With the manager
with gw.config.update(ref_image=l8_224077_20200518_B2):
    with gw.open(l8_224078_20200518) as src:
        print(src.gw.bounds)
(717345.0, -2832795.0, 778575.0, -2776995.0)
(694005.0, -2812065.0, 754185.0, -2766615.0)
(694005.0, -2812065.0, 754185.0, -2766615.0)

Reference settings: Sensors#

Because rasters are opened as DataArrays, the band coordinates will be named. By default, the bands will be named by their index position (starting at 1). It might, however, be more intuitive to store the band names as strings, where the names correspond to the sensor wavelengths. In GeoWombat, you can set the band names explicitly upon opening a file by using the :func:geowombat.open band_names keyword. Alternatively, if the sensor is known (and supported by GeoWombat), then you can set the band names by specifying the sensor name in the configuration settings.

Note

In the example below, the example raster comes from a Landsat image. However, only the visible (blue, green, and red) wavelengths are stored. Thus, we use ‘rgb’ as the sensor name. If we had a full 6-band Landsat 7 image, for example, we could use the ‘l7’ sensor flag.

import geowombat as gw
from geowombat.data import l8_224078_20200518

# Without the manager
with gw.open(l8_224078_20200518) as src:
    print(src.band)

# With the manager
with gw.config.update(sensor='bgr'):
    with gw.open(l8_224078_20200518) as src:
        print(src.band)
<xarray.DataArray 'band' (band: 3)>
array([1, 2, 3])
Coordinates:
  * band     (band) int64 1 2 3
<xarray.DataArray 'band' (band: 3)>
array(['blue', 'green', 'red'], dtype='<U5')
Coordinates:
  * band     (band) <U5 'blue' 'green' 'red'

To see all available sensor names, use the avail_sensors property.

with gw.open(l8_224078_20200518) as src:
    for sensor_name in src.gw.avail_sensors:
        print(sensor_name)

For a short description of the sensor, use the sensor_names property.

with gw.open(l8_224078_20200518) as src:
    for sensor_name, description in src.gw.sensor_names.items():
        print('{}: {}'.format(sensor_name.ljust(15), description))
rgb            : red, green, and blue
rgbn           : red, green, blue, and NIR
bgr            : blue, green, and red
bgrn           : blue, green, red, and NIR
l5             : Landsat 5 Thematic Mapper (TM)
l7             : Landsat 7 Enhanced Thematic Mapper Plus (ETM+) without panchromatic and thermal bands
l7th           : Landsat 7 Enhanced Thematic Mapper Plus (ETM+) with thermal band
l7mspan        : Landsat 7 Enhanced Thematic Mapper Plus (ETM+) with panchromatic band
l7pan          : Landsat 7 panchromatic band
l8             : Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) without panchromatic and thermal bands
l8l7           : Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with 6 Landsat 7-like bands
l8l7mspan      : Landsat 8 Operational Land Imager (OLI) and panchromatic band with 6 Landsat 7-like bands
l8th           : Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with thermal band
l8pan          : Landsat 8 panchromatic band
l9             : Landsat 9 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) without panchromatic and thermal bands
l9l7           : Landsat 9 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with 6 Landsat 7-like bands
l9l7mspan      : Landsat 9 Operational Land Imager (OLI) and panchromatic band with 6 Landsat 7-like bands
l9th           : Landsat 9 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with thermal band
l9pan          : Landsat 9 panchromatic band
s2             : Sentinel 2 Multi-Spectral Instrument (MSI) without 3 60m bands (coastal, water vapor, cirrus)
s2a            : Sentinel 2 Multi-Spectral Instrument (MSI) without 3 60m bands (coastal, water vapor, cirrus)
s2b            : Sentinel 2 Multi-Spectral Instrument (MSI) without 3 60m bands (coastal, water vapor, cirrus)
s2f            : Sentinel 2 Multi-Spectral Instrument (MSI) with 3 60m bands (coastal, water vapor, cirrus)
s2l7           : Sentinel 2 Multi-Spectral Instrument (MSI) with 6 Landsat 7-like bands
s2al7          : Sentinel 2 Multi-Spectral Instrument (MSI) with 6 Landsat 7-like bands
s2bl7          : Sentinel 2 Multi-Spectral Instrument (MSI) with 6 Landsat 7-like bands
s210           : Sentinel 2 Multi-Spectral Instrument (MSI) with 4 10m (visible + NIR) bands
s220           : Sentinel 2 Multi-Spectral Instrument (MSI) with 6 20m bands
s2cloudless    : Sentinel 2 Multi-Spectral Instrument (MSI) with 10 bands for s2cloudless
ps             : PlanetScope with 4 (visible + NIR) bands
qb             : Quickbird with 4 (visible + NIR) bands
ik             : IKONOS with 4 (visible + NIR) bands
mcd43a4        : MODIS Nadir BRDF-Adjusted Reflectance Daily 500m with 7 bands

The following is a list of all available sensor names. This documentation may become out of date, if so please refer to geowombat/core/properties.py for the full list.

Abbreviated Name

Description

‘rgb’

red, green, and blue

‘rgbn’

red, green, blue, and NIR

‘bgr’

blue, green, and red

‘bgrn’

blue, green, red, and NIR

‘l5’

Landsat 5 Thematic Mapper (TM)

‘l7’

Landsat 7 Enhanced Thematic Mapper Plus (ETM+) without panchromatic and thermal bands

‘l7th’

Landsat 7 Enhanced Thematic Mapper Plus (ETM+) with thermal band

‘l7mspan’

Landsat 7 Enhanced Thematic Mapper Plus (ETM+) with panchromatic band

‘l7pan’

Landsat 7 panchromatic band

‘l8’

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) without panchromatic and thermal bands

‘l8l7’

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with 6 Landsat 7-like bands

‘l8l7mspan’

Landsat 8 Operational Land Imager (OLI) and panchromatic band with 6 Landsat 7-like bands

‘l8th’

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with thermal band

‘l8pan’

Landsat 8 panchromatic band

‘s2’

Sentinel 2 Multi-Spectral Instrument (MSI) without 3 60m bands (coastal, water vapor, cirrus)

‘s2f’

Sentinel 2 Multi-Spectral Instrument (MSI) with 3 60m bands (coastal, water vapor, cirrus)

‘s2l7’

Sentinel 2 Multi-Spectral Instrument (MSI) with 6 Landsat 7-like bands

‘s210’

Sentinel 2 Multi-Spectral Instrument (MSI) with 4 10m (visible + NIR) bands

‘s220’

Sentinel 2 Multi-Spectral Instrument (MSI) with 6 20m bands

‘s2cloudless’

Sentinel 2 Multi-Spectral Instrument (MSI) with 10 bands for s2cloudless

‘ps’

PlanetScope with 4 (visible + NIR) bands

‘qb’

Quickbird with 4 (visible + NIR) bands

‘ik’

IKONOS with 4 (visible + NIR) bands