geowatch.utils.util_gdal module

SeeAlso

util_raster.py

References

https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-multi https://gis.stackexchange.com/a/241810 https://trac.osgeo.org/gdal/wiki/UserDocs/GdalWarp#WillincreasingRAMincreasethespeedofgdalwarp https://github.com/OpenDroneMap/ODM/issues/778

Todo

TODO test this and see if it’s safe to add:

–config GDAL_PAM_ENABLED NO

Removes .aux.xml sidecar files and puts them in the geotiff metadata ex. histogram from fmask https://stackoverflow.com/a/51075774 https://trac.osgeo.org/gdal/wiki/ConfigOptions#GDAL_PAM_ENABLED https://gdal.org/drivers/raster/gtiff.html#georeferencing

geowatch.utils.util_gdal.import_gdal()[source]
geowatch.utils.util_gdal.import_osr()[source]
class geowatch.utils.util_gdal.DummyLogger[source]

Bases: object

warning(msg, *args)[source]
class geowatch.utils.util_gdal.GDalCommandBuilder(command_name)[source]

Bases: object

Helper for templating out gdal commands.

Usage is to give the builder the name of the main command, and then you can add flags, single-value options, multi-value options, and positional arguments. Convinience methods exist for common option profiles. When done call finalize to get the final command string.

Example

>>> from geowatch.utils.util_gdal import *  # NOQA
>>> builder = GDalCommandBuilder('madeupcmd')
>>> builder.flags.append('-flag1')
>>> builder.flags.append('-flag2')
>>> builder['-of'] = 'VRT'
>>> builder['-f1'] = 'flag1'
>>> builder['-f2'] = 'a b c d'
>>> builder['-xywh'] = ['x', 'y', 'w', 'h']
>>> builder.options['-oo']['FOO'] = 'bar'
>>> builder.options['-oo']['BAZ'] = 'biz'
>>> builder.options['--config']['CPL_LOG'] = '/dev/null'
>>> builder.options['--config']['SPACE_PATH'] = '/path with/spaces'
>>> builder.append('pos1')
>>> builder.append('pos2')
>>> builder.set_cog_options()
>>> command = builder.finalize()
>>> print(command)
property config_options
property common_options
property options
property flags
append(item)[source]
extend(items)[source]
finalize()[source]
Returns:

the final bash command

Return type:

str

set_cog_options(compress='DEFLATE', blocksize=256, overviews='AUTO', overview_resampling='CUBIC')[source]
geowatch.utils.util_gdal.gdal_single_translate(in_fpath, out_fpath, pixel_box=None, blocksize=256, compress='DEFLATE', tries=1, cooldown=1, backoff=1.0, verbose=0, eager=True, gdal_cachemax=None, num_threads=None, use_tempfile=True, timeout=None, overviews='AUTO', overview_resampling='CUBIC')[source]

Crops geotiffs using pixels

Parameters:

  • in_fpath (PathLike) – geotiff to translate

  • out_fpath (PathLike) – output geotiff

  • pixel_box (kwimage.Boxes) – box to crop to in pixel space.

  • blocksize (int) – COG tile size

  • compress (str) – gdal compression

  • verbose (int) – verbosity level

  • cooldown (int) – number of seconds to wait (i.e. “cool-down”) between tries

  • backoff (float) – multiplier applied to cooldown between attempts. default: 1 (no backoff)

  • eager (bool) – if True, executes the command, if False returns a list of all the bash commands that would be executed, suitable for use in a command queue.

  • gdal_cachemax (int | None) – i/o block cache size in megabytes.

  • timeout (None | float) – number of seconds allowed to run gdal before giving up

CommandLine

xdoctest -m geowatch.utils.util_gdal gdal_single_translate

Example

>>> from geowatch.utils.util_gdal import *  # NOQA
>>> from geowatch.utils.util_gdal import _demo_geoimg_with_nodata
>>> from geowatch.utils import util_gis
>>> from geowatch.gis import geotiff
>>> in_fpath = ub.Path(_demo_geoimg_with_nodata())
>>> info = geotiff.geotiff_crs_info(in_fpath)
>>> gdf = util_gis.crs_geojson_to_gdf(info['geos_corners'])
>>> gdf = gdf.to_crs(4326) # not really wgs84, this is crs84
>>> # Test CRS84 cropping
>>> crs84_poly = kwimage.Polygon.coerce(gdf['geometry'].iloc[0])
>>> crs84_epsg = gdf.crs.to_epsg()
>>> crs84_space_box = crs84_poly.scale(0.5, about='center').to_boxes()
>>> crs84_out_fpath = ub.augpath(in_fpath, prefix='_crs84_crop')
>>> gdal_single_warp(in_fpath, crs84_out_fpath, local_epsg=crs84_epsg, space_box=crs84_space_box)

>>> # Test UTM cropping
>>> utm_gdf = util_gis.project_gdf_to_local_utm(gdf)
>>> utm_poly = kwimage.Polygon.coerce(utm_gdf['geometry'].iloc[0])
>>> utm_epsg = utm_gdf.crs.to_epsg()
>>> utm_space_box = utm_poly.scale(0.5, about='center').to_boxes()
>>> utm_out_fpath = ub.augpath(in_fpath, prefix='_utmcrop')
>>> gdal_single_warp(in_fpath, utm_out_fpath, local_epsg=utm_epsg, space_box=utm_space_box, box_epsg=utm_epsg)

>>> # Test Pixel cropping
>>> pxl_poly = kwimage.Polygon(exterior=info['pxl_corners'])
>>> pixel_box = pxl_poly.scale(0.5, about='center').to_boxes()
>>> pxl_out_fpath = ub.augpath(in_fpath, prefix='_pxlcrop')
>>> gdal_single_translate(in_fpath, pxl_out_fpath, pixel_box=pixel_box)

>>> # xdoctest: +REQUIRES(--show)
>>> import kwplot
>>> kwplot.autompl()
>>> imdata0 = kwimage.normalize(kwimage.imread(in_fpath, nodata='float'))
>>> imdata1 = kwimage.normalize(kwimage.imread(crs84_out_fpath, nodata='float'))
>>> imdata2 = kwimage.normalize(kwimage.imread(utm_out_fpath, nodata='float'))
>>> imdata3 = kwimage.normalize(kwimage.imread(pxl_out_fpath, nodata='float'))
>>> kwplot.imshow(imdata0, pnum=(1, 4, 1), title='orig')
>>> kwplot.imshow(imdata1, pnum=(1, 4, 2), title='crs84-crop')
>>> kwplot.imshow(imdata2, pnum=(1, 4, 3), title='utm-crop')
>>> kwplot.imshow(imdata3, pnum=(1, 4, 4), title='pxl-crop')

Example

>>> from geowatch.utils.util_gdal import *  # NOQA
>>> import kwimage
>>> in_fpath = kwimage.grab_test_image_fpath('amazon')
>>> out_fpath = 'foo.tif'
>>> commands = gdal_single_translate(in_fpath, out_fpath, eager=False)
>>> assert len(commands) == 2
geowatch.utils.util_gdal.gdal_single_warp(in_fpath, out_fpath, space_box=None, local_epsg=4326, box_epsg=4326, nodata=None, rpcs=None, blocksize=256, compress='DEFLATE', overviews='AUTO', overview_resampling='CUBIC', as_vrt=False, use_te_geoidgrid=False, dem_fpath=None, error_logfile=None, tries=1, verbose=0, force_spatial_res=None, eager=True, timeout=None, cooldown=1, backoff=1.0, use_tempfile=True, gdal_cachemax=None, num_threads=None, warp_memory=None)[source]

Wrapper around gdalwarp

Parameters:

  • in_fpath (PathLike) – input geotiff path

  • out_fpath (PathLike) – output geotiff path

  • space_box (kwimage.Boxes) – Should be traditional crs84 ltrb (or lbrt?) – i.e. (lonmin, latmin, lonmax, latmax) - when box_epsg is 4326

  • local_epsg (int) – EPSG code for the CRS the final geotiff will be projected into. This should be the UTM zone for the region if known. Otherwise It can be 4326 to project into WGS84 or CRS84 (not sure which axis ordering it will use by default).

  • box_epsg (int) – this is the EPSG of the bounding box. Should usually be 4326.

  • nodata (int | None) – only specify if in_fpath does not already have a nodata value

  • rpcs (dict) – the “rpc_transform” from geowatch.gis.geotiff.geotiff_crs_info, if that information is available and orthorectification is desired.

  • as_vrt (bool) – undocumented

  • use_te_geoidgrid (bool) – undocumented

  • dem_fpath (bool) – undocumented

  • error_logfile (None | PathLike) – If specified, errors will be logged to this filepath.

  • tries (int) – gdal can be flakey, set to force some number of retries. Ignored if eager is False.

  • force_spatial_res (float | tuple(float, float)) – Force spatial resolution for output images.

  • eager (bool) – if True, executes the command, if False returns a list of all the bash commands that would be executed, suitable for use in a command queue.

  • cooldown (int) – number of seconds to wait (i.e. “cool-down”) between tries

  • backoff (float) – multiplier applied to cooldown between attempts. default: 1 (no backoff)

  • timeout (None | float) – number of seconds allowed to run gdal before giving up

Returns:

Nothing if executing the command. If eager=False, returns the commands that would have been executed.

Return type:

None | List[str]

Notes

In gdalwarp:

-s_srs - Set source spatial reference -t_srs - Set target spatial reference

-te_srs - Specifies the SRS in which to interpret the coordinates given with -te. -te - Set georeferenced extents of output file to be created

Example

>>> import kwimage
>>> from geowatch.utils.util_gdal import gdal_single_warp
>>> # Note: this is a network test that depends on external resources
>>> # we may want to disable, or try to make more robust.
>>> in_fpath = '/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/23/K/PQ/2019/6/S2B_23KPQ_20190623_0_L2A/B03.tif'
>>> from osgeo import gdal
>>> info = gdal.Info(in_fpath, format='json')
>>> bound_poly = kwimage.Polygon.coerce(info['wgs84Extent'])
>>> crop_poly = bound_poly.scale(0.02, about='centroid')
>>> space_box = crop_poly.to_boxes()
>>> out_fpath = ub.Path.appdir('geowatch/doctests/util_gdal').ensuredir() / 'cropped.tif'
>>> error_logfile = '/dev/null'
>>> gdal_single_warp(in_fpath, out_fpath, space_box=space_box, error_logfile=error_logfile, verbose=3)
>>> # xdoctest: +REQUIRES(--show)
>>> import kwplot
>>> kwplot.autompl()
>>> data = kwimage.imread(out_fpath)
>>> canvas = kwimage.normalize_intensity(data)
>>> kwplot.imshow(canvas)

Example

>>> # Test non-eager version
>>> import kwimage
>>> from geowatch.utils.util_gdal import gdal_single_warp
>>> in_fpath = '/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/23/K/PQ/2019/6/S2B_23KPQ_20190623_0_L2A/B03.tif'
>>> from osgeo import gdal
>>> info = gdal.Info(in_fpath, format='json')
>>> bound_poly = kwimage.Polygon.coerce(info['wgs84Extent'])
>>> crop_poly = bound_poly.scale(0.02, about='centroid')
>>> space_box = crop_poly.to_boxes()
>>> out_fpath = ub.Path.appdir('geowatch/doctests/util_gdal').ensuredir() / 'cropped.tif'
>>> commands = gdal_single_warp(in_fpath, out_fpath, space_box=space_box, verbose=3, eager=False)
>>> assert len(commands) == 2

References

https://gdal.org/programs/gdalwarp.html

geowatch.utils.util_gdal.gdal_multi_warp(in_fpaths, out_fpath, space_box=None, local_epsg=4326, box_epsg=4326, nodata=None, tries=1, cooldown=1, backoff=1.0, blocksize=256, compress='DEFLATE', overviews='AUTO', overview_resampling='CUBIC', error_logfile=None, _intermediate_vrt=False, verbose=0, return_intermediate=False, force_spatial_res=None, eager=True, gdal_cachemax=None, num_threads=None, warp_memory=None, use_tempfile=True, timeout=None, **kwargs)[source]

See gdal_single_warp() for args

NOTE: it is important to set the nodata argument for gdalmerge [SO187522] if you want to preserver them.

Example

>>> from geowatch.utils.util_gdal import *  # NOQA
>>> from geowatch.utils.util_gdal import _demo_geoimg_with_nodata
>>> from geowatch.utils import util_gis
>>> from geowatch.gis import geotiff
>>> in_fpath = ub.Path(_demo_geoimg_with_nodata())
>>> info = geotiff.geotiff_crs_info(in_fpath)
>>> gdf = util_gis.crs_geojson_to_gdf(info['geos_corners'])
>>> gdf = gdf.to_crs(4326) # not really wgs84, this is crs84
>>> crs84_poly = kwimage.Polygon.coerce(gdf['geometry'].iloc[0])
>>> crs84_epsg = gdf.crs.to_epsg()
>>> crs84_space_box = crs84_poly.scale(0.5, about='center').to_boxes()
>>> crs84_out_fpath = ub.augpath(in_fpath, prefix='_crs84_crop_2_')
>>> in_fpaths = [in_fpath, in_fpath, in_fpath]
>>> # Test that timeouts work
>>> import pytest
>>> with pytest.raises(subprocess.TimeoutExpired):
...     gdal_multi_warp(in_fpaths, crs84_out_fpath, local_epsg=crs84_epsg, space_box=crs84_space_box, timeout=0)

# Dev case: # The first part should work, but later parts should fail if 0:

with ub.Timer() as timer:

gdal_multi_warp(in_fpaths, crs84_out_fpath, local_epsg=crs84_epsg, space_box=crs84_space_box, verbose=3)

print(f’timer.elapsed={timer.elapsed}’) gdal_multi_warp(in_fpaths, crs84_out_fpath, local_epsg=crs84_epsg, space_box=crs84_space_box, timeout=timer.elapsed * 0.6, verbose=3)

Example

>>> # xdoctest: +REQUIRES(--slow)
>>> # Uses data from the data cube with extra=1
>>> from geowatch.utils.util_gdal import *  # NOQA
>>> import ubelt as ub
>>> local_epsg = 32635
>>> space_box = kwimage.Polygon.random().bounding_box().to_ltrb()
>>> out_fpath = 'lazy_multi_warp.tif'
>>> in_fpaths = ['dummy1.tif', 'dummy2.tif']
>>> commands = gdal_multi_warp(
>>>     in_fpaths, out_fpath=out_fpath, space_box=space_box,
>>>     local_epsg=local_epsg, verbose=3, eager=False)
>>> for cmd in commands:
>>>     print(cmd)
Returns:

Nothing if executing the command. If eager=False, returns the commands that would have been executed.

Return type:

None | List[str]

References

geowatch.utils.util_gdal.list_gdal_drivers()[source]

List all drivers currently available to GDAL to create a raster

Returns:

list((driver_shortname, driver_longname, list(driver_file_extension)))

Example

>>> from geowatch.utils.util_gdal import *
>>> drivers = list_gdal_drivers()
>>> print('drivers = {}'.format(ub.urepr(drivers, nl=1)))
>>> assert ('GTiff', 'GeoTIFF', ['tif', 'tiff']) in drivers
geowatch.utils.util_gdal.GdalOpen(path, mode='r', **kwargs)[source]

A simple context manager for friendlier gdal use.

Returns:

GdalDataset

Example

>>> # xdoctest: +REQUIRES(--network)
>>> from geowatch.utils.util_gdal import *
>>> from osgeo import gdal
>>> from geowatch.demo.landsat_demodata import grab_landsat_product
>>> path = grab_landsat_product()['bands'][0]
>>> #
>>> # standard use:
>>> dataset = gdal.Open(path)
>>> print(dataset.GetDescription())  # do stuff
>>> del dataset  # or 'dataset = None'
>>> #
>>> # equivalent:
>>> with GdalOpen(path) as dataset:
>>>     print(dataset.GetDescription())  # do stuff
>>> #
>>> # open for writing:
>>> with GdalOpen(path, gdal.GA_Update) as dataset:
>>>     print(dataset.GetDescription())  # do stuff
class geowatch.utils.util_gdal.GdalSupressWarnings[source]

Bases: object

References

https://gdal.org/api/python_gotchas.html

This currently will suppress warnings entirely. We could build this into something a big nicer.

handler(err_level, err_no, err_msg)[source]
class geowatch.utils.util_gdal.GdalDataset(_GdalDataset__ref, _path='?', _str_mode='?')[source]

Bases: NiceRepr

A wrapper around gdal.Open and the underlying dataset it returns.

This object is completely transparent and offers the same API as the osgeo.gdal.Dataset returned by :func`:osgeo.gdal.GDalOpen`.

This object can be used as a context manager. By default the GDAL dataset is opened when the object is created, and it is closed when either close is called or the __exit__ method is called by a context manager. When the object is closed the underlying GDAL objet is dereferenced and garbage collected.

Parameters:

  • path (PathLike) – a path or string referencing a gdal image file

  • mode (str | int) – a gdal GA (Gdal Access) integer code or a string that can be: ‘readonly’ or ‘update’ or the equivalent standard mode codes: ‘r’ and ‘w+’.

  • virtual_retries (int) – If the path is a reference to a virtual file system (i.e. starts with vsi) then we try to open it this many times before we finally fail.

Example

>>> # Demonstrate use cases of this object
>>> from geowatch.utils.util_gdal import *
>>> import kwimage
>>> # Grab demo path we can test with
>>> path = kwimage.grab_test_image_fpath()
>>> #
>>> #
>>> # Method1: Use GDalOpen exactly the same as gdal.Open
>>> ref = GdalDataset.open(path)
>>> print(f'{ref=!s}')
>>> assert not ref.closed
>>> ref.GetDescription()  # use GDAL API exactly as-is
>>> assert not ref.closed
>>> ref.close()  # Except you can now do this
>>> print(f'{ref=!s}')
>>> assert ref.closed
>>> #
>>> #
>>> # Method2: Use GDalOpen exactly the same as gdal.GdalDataset
>>> with GdalDataset.open(path, mode='r') as ref:
>>>     ref.GetDescription()  # do stuff
>>>     print(f'{ref=!s}')
>>>     assert not ref.closed
>>> print(f'{ref=!s}')
>>> assert ref.closed

Example

>>> # Test virtual filesystem
>>> # xdoctest: +REQUIRES(--network)
>>> from geowatch.utils.util_gdal import *  # NOQA
>>> path = '/vsicurl/https://i.imgur.com/KXhKM72.png'
>>> ref = GdalDataset.open(path)
>>> data = ref.GetRasterBand(1).ReadAsArray()
>>> assert data.sum() == 37109758

Do not call this method directly. Use GdalDataset.open

classmethod open(path, mode='r', virtual_retries=3, cooldown=0.1)[source]

Create a new dataset in read or write mode.

Parameters:

  • path (str | PathLike) – the path or URI to open

  • mode (str) – either ‘r’ for read only or ‘w+’ for update mode.

  • virtual_retries (int) – number of times to attempt to open the dataset when the URI points to a non-local resource

  • cooldown (float) – seconds between retry attempts

Returns:

GdalDataset

classmethod coerce(data, mode=None, **kwargs)[source]

Ensures the underlying object is a gdal dataset.

property closed
property mode
property shape
close()[source]

Closes this dataset.

Part of the GDalOpen Wrapper. Closes this dataset and dereferences the underlying GDAL object.

Note: this will not work if the __ref attribute as accessed outside of this wrapper class.

info()[source]

Information similar to gdalinfo (in json format)

get_overview_info()[source]

get number of overviews for each band.

get_all_metadata()[source]

References

https://gdal.org/user/raster_data_model.html