geowatch.gis.geotiff module¶

Tools to work with geotiff metadata.

geowatch.gis.geotiff.geotiff_metadata(gpath, elevation='gtop30', strict=False, supress_warnings=False)[source]¶

Extract all relevant metadata we know how to extract.

Parameters:

gpath (str) – path to the geotiff of interest
elevation (str) – method for extracting elevation data.

Example

>>> # xdoctest: +REQUIRES(env:SLOW_DOCTEST)
>>> # xdoctest: +REQUIRES(--network)
>>> from geowatch.gis.geotiff import *  # NOQA
>>> url = ('http://storage.googleapis.com/gcp-public-data-landsat/'
...        'LC08/01/044/034/LC08_L1GT_044034_20130330_20170310_01_T2/'
...        'LC08_L1GT_044034_20130330_20170310_01_T2_B11.TIF')
>>> gpath = ub.grabdata(url, appname='geowatch')
>>> info = geotiff_metadata(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))

>>> import ubelt as ub
>>> url = ('http://storage.googleapis.com/gcp-public-data-landsat/'
...        'LC08/01/037/029/LC08_L1TP_037029_20130602_20170310_01_T1/'
...        'LC08_L1TP_037029_20130602_20170310_01_T1_B2.TIF')
>>> gpath = ub.grabdata(url, appname='geowatch')

>>> info = geotiff_metadata(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))

geowatch.gis.geotiff.geotiff_header_info(gpath_or_ref)[source]¶

Extract relevant metadata information from a geotiff header.

Example

>>> # xdoctest: +REQUIRES(env:SLOW_DOCTEST)
>>> from geowatch.gis.geotiff import *  # NOQA
>>> from geowatch.demo.dummy_demodata import dummy_rpc_geotiff_fpath
>>> from geowatch.demo.landsat_demodata import grab_landsat_product
>>> gpath_or_ref = gpath = dummy_rpc_geotiff_fpath()
>>> info = geotiff_header_info(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))
>>> gpath_or_ref = gpath = grab_landsat_product()['bands'][0]
>>> info = geotiff_header_info(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))

geowatch.gis.geotiff.geotiff_crs_info(gpath_or_ref, force_affine=False, elevation='gtop30', verbose=0)[source]¶

Use GDAL to extract coordinate system information about a geo_tiff.

Builds transformations between pixel, geotiff-world, utm, and wgs84 spaces

Parameters:

gpath (str) – path to the image file
force_affine (bool) – if True ignores RPC information
elevation (str) – method used to determine the elevation when RPC information is used. Available options are:
- “open-elevation”
- “gtop30”

Example

>>> # xdoctest: +REQUIRES(env:SLOW_DOCTEST)
>>> from geowatch.gis.geotiff import *  # NOQA
>>> from geowatch.demo.dummy_demodata import dummy_rpc_geotiff_fpath
>>> gpath_or_ref = gpath = dummy_rpc_geotiff_fpath()
>>> info = geotiff_crs_info(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1, sort=False)))
>>> assert info['is_rpc']
>>> assert info['img_shape'] == (2000, 2000)

>>> # xdoctest: +REQUIRES(--network)
>>> gpath_or_ref = gpath = ub.grabdata(
>>>     'https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif',
>>>     appname='geowatch/demodata', hash_prefix='10a2ebcdcd95582')
>>> info = geotiff_crs_info(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1, sort=False)))
>>> assert not info['is_rpc']
>>> assert info['img_shape'] == (515, 514)

>>> # The public gateways seem to be too slow in serving the content
>>> # xdoctest: +REQUIRES(--ipfs)
>>> from geowatch.demo.nitf_demodata import grab_nitf_fpath
>>> gpath_or_ref = gpath = grab_nitf_fpath('i_3004g.ntf')
>>> info = geotiff_crs_info(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1, sort=False)))
>>> assert not info['is_rpc']
>>> assert info['img_shape'] == (512, 512)

tf = info[‘wgs84_to_wld’]

geowatch.gis.geotiff.make_crs_info_object(osr_crs)[source]¶

Parameters:: osr_crs (osr.SpatialReference) – an osr object from gdal

Example

>>> from geowatch.gis.geotiff import *  # NOQA
>>> from geowatch.utils import util_gdal
>>> osr = util_gdal.import_osr()
>>> osr_crs = osr.SpatialReference()
>>> osr_crs.ImportFromEPSG(4326)
>>> crs_info = make_crs_info_object(osr_crs)
>>> print('crs_info = {}'.format(ub.urepr(crs_info, nl=1)))
>>> osr_crs = osr.SpatialReference()
>>> osr_crs.ImportFromEPSG(4326)
>>> osr_crs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
>>> crs_info = make_crs_info_object(osr_crs)
>>> print('crs_info = {}'.format(ub.urepr(crs_info, nl=1)))
>>> osr_crs.ImportFromEPSG(32744)
>>> crs_info = make_crs_info_object(osr_crs)
>>> print('crs_info = {}'.format(ub.urepr(crs_info, nl=1)))

geowatch.gis.geotiff.axis_mapping_int_to_text(axis_mapping_int)[source]¶

References

https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order

Notes

OAMS_TRADITIONAL_GIS_ORDER means that for geographic CRS with
lat/long order, the data will still be long/lat ordered. Similarly for a projected CRS with northing/easting order, the data will still be easting/northing ordered
OAMS_AUTHORITY_COMPLIANT means that the data axis will be
identical to the CRS axis. This is the default value when instantiating OGRSpatialReference
OAMS_CUSTOM means that the data axis are customly defined with
SetDataAxisToSRSAxisMapping

geowatch.gis.geotiff.memo_auth_from_wkt(wkt)[source]¶: This benchmarks as an expensive operation, memoize it.

exception geowatch.gis.geotiff.InvalidFormat[source]¶: Bases: Exception

geowatch.gis.geotiff.geotiff_filepath_info(gpath, fast=True)[source]¶

Attempt to parse information out of a path to a geotiff file.

Information provided here is purely heuristic. Generally filepath information is not robust.

Several huerstics are currently implemented for:

Sentinel 2
Landsat-8
WorldView3

See [S2_Name_2016] and [S3_Name].

Parameters:

gpath (str) – a path to an image that uses a standard naming convention (may include subdirectories that contain relevant information) .
fast (bool) – if True stops when a hueristic matches well enough, otherwise tries multiple hueristics. Defaults to True.

SeeAlso:

parse_landsat_product_id - specific to the landsat spec

Example

>>> from geowatch.gis.geotiff import *  # NOQA
>>> inputs = [
>>>     '18JUN18071532-S3DMR1C2.NTF',
>>>     'LC08_L1GT_029030_20151209_20160131_01_RT',
>>>     'LC08_L1TP_162043_20130401_20170505_01_T1_B8.tif',
>>>     'LC08_L1GT_017037_20190303_20190309_01_T2_B8.TIF',
>>>     '03AUG16WV031000016AUG03125811-P1BS-500827464010_01_P002_________AAE_0AAAAABPABB0.NTF',
>>>     '13APR16WV030900016APR13131153-P1BS_R1C1-500659828010_01_P001____AAE_0AAAAABPABM0.NTF.NTF',
>>>     '011777481010_01_P001_MUL/17SEP07021826-M1BS-011777481010_01_P001.TIF',
>>>     '011778213010_01_P001_PAN/14FEB03022210-P1BS-011778213010_01_P001.NTF',
>>>     'T17SMS_20151020T162042_TCI.jp2',
>>>     'T39QXG_20160711T065622_TCI.jp2',
>>>     'S2B_MSIL1C_20181219T022109_N0207_R003_T52SDG_20181219T051028.SAFE/GRANULE/L1C_T52SDG_A009324_20181219T022640/IMG_DATA/T52SDG_20181219T022109_TCI.jp2',
>>>     'S2B_MSIL1C_20181219T022109_N0207_R003_T52SDG_20181219T051028/S2B_MSIL1C_20181219T022109_N0207_R003_T52SDG_20181219T051028.SAFE/GRANULE/L1C_T52SDG_A009324_20181219T022640/IMG_DATA/T52SDG_20181219T022109_B01.jp2',
>>>     'S2A_MSIL1C_20151021T022702_N0204_R003_T52SDG_20151021T022701.SAFE/GRANULE/L1C_T52SDG_A001716_20151021T022701/IMG_DATA/T52SDG_20151021T022702_TCI.jp2',
>>>     'LC08_L2SP_217076_20190107_20211008_02_T1_T23KPQ_B1.tif',
>>>     'S2B_MSI_L2A_T23KPQ_20190114_20211008_SR_B01.tif',
>>>     'S2A_MSI_L2A_T39RVK_20180803_20211102_SR_SOZ4.tif',
>>> ]
>>> gpath = inputs[-1]
>>> print('gpath = {!r}'.format(gpath))
>>> for gpath in inputs:
>>>     print('----')
>>>     info = geotiff_filepath_info(gpath)
>>>     print('gpath = {}'.format(ub.urepr(gpath, nl=1)))
>>>     print('info = {}'.format(ub.urepr(info, nl=2)))
>>>     if len(info['sensor_candidates']) == 0:
>>>         print(ub.color_text('NO HUERISTIC', 'red'))
>>>     else:
>>>         assert len(info['sensor_candidates']) == len(set(info['sensor_candidates']))
>>>         assert info['filename_meta']

Example

>>> from geowatch.gis.geotiff import *  # NOQA
>>> gpath = 'LC08_L1TP_037029_20130602_20170310_01_T1_B1.TIF'
>>> info = geotiff_filepath_info(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))
>>> assert info['filename_meta']['sensor_code'] == 'C'
>>> assert info['filename_meta']['sat_code'] == '08'
>>> assert info['filename_meta']['sat_code'] == '08'
>>> assert info['filename_meta']['collection_category'] == 'T1'
>>> assert info['filename_meta']['suffix'] == 'B1'

Example

>>> # xdoctest: +REQUIRES(--network)
>>> # Test extact info from real landsat product files
>>> from geowatch.demo.landsat_demodata import grab_landsat_product  # NOQA
>>> from geowatch.gis.geotiff import *  # NOQA
>>> product = grab_landsat_product()
>>> band_infos = [geotiff_filepath_info(gpath) for gpath in product['bands']]
>>> meta_infos = [geotiff_filepath_info(gpath) for gpath in product['meta'].values()]
>>> assert not any(d is None for d in band_infos)
>>> assert not any(d is None for d in meta_infos)

References

[S2_Name_2016]

https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/naming-convention

[S3_Name]

https://sentinel.esa.int/web/sentinel/user-guides/sentinel-3-altimetry/naming-conventions

geowatch.gis.geotiff.parse_sentinel2_product_id(parts)[source]¶

Try to parse the Sentinel-2 pre-2016 and post-2016 safedir formats.

Note that unlike parse_landsat_product_id, which expects a band file basename, this presently purports to plurally parse pieces of path postfixedly (it parses the whole path, backwards :))

TODO extend this to parsing the granuledir and band file formats as well. For now, we just need all names to be minimally parseable, even if some info is incorrect.

General plan is to check the old formats strictly first, and then check the new safedir loosely as a default

Example

parts = [‘S2B_MSI_L2A_T23KPQ_20190114_20211008_SR_B01.tif’] parts = [‘S2A_MSI_L2A_T39RVK_20180803_20211102_SR_SOZ4.tif’] parse_sentinel2_product_id(parts)

geowatch.gis.geotiff.parse_landsat_product_id(product_id)[source]¶

Extract information from a landsat produt id

See [LanSatName], [LS_578], [ExampleLandSat], [LandSatSuffixFormat], [LandsatProcLevels], [LandsatL2Names], and [LandSatARDDocs].

Parameters:: product_id (str) – this is typically the filename (without extension!) of a landsat product, as described in [LanSatName].

Example

>>> from geowatch.gis.geotiff import *  # NOQA
>>> product_id = 'LC08_L1TP_037029_20130602_20170310_01_T1'
>>> ls_meta = parse_landsat_product_id(product_id)

Example

>>> from geowatch.gis.geotiff import *  # NOQA
>>> gpath = 'LC08_L1TP_037029_20130602_20170310_01_T1_B1'
>>> info = parse_landsat_product_id(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))
>>> assert info['sensor_code'] == 'C'
>>> assert info['sat_code'] == '08'
>>> assert info['sat_code'] == '08'
>>> assert info['collection_category'] == 'T1'
>>> assert info['suffix'] == 'B1'
>>> assert info['band_num'] == 1

>>> from geowatch.gis.geotiff import *  # NOQA
>>> gpath = 'LC08_CU_029005_20181208_20210503_02_QA_LINEAGE.TIF'
>>> info = parse_landsat_product_id(gpath)
>>> print('info = {}'.format(ub.urepr(info, nl=1)))

Example

>>> # xdoctest: +REQUIRES(--network)
>>> # Test on real landsat data
>>> from geowatch.demo.landsat_demodata import grab_landsat_product  # NOQA
>>> from geowatch.gis.geotiff import *  # NOQA
>>> product = grab_landsat_product()
>>> band_prodids = [ub.augpath(gpath, dpath='', ext='') for gpath in product['bands']]
>>> band_infos = [parse_landsat_product_id(product_id) for product_id in band_prodids]
>>> assert ub.allsame([ub.dict_diff(d, ['band_num', 'suffix', 'channels', 'band']) for d in band_infos])
>>> meta_prodids = [ub.augpath(gpath, dpath='', ext='') for gpath in product['meta'].values()]
>>> meta_infos = [parse_landsat_product_id(product_id) for product_id in meta_prodids]
>>> assert ub.allsame([ub.dict_diff(d, ['suffix', 'channels', 'band']) for d in meta_infos])

References

[LanSatName] (1,2)

https://www.usgs.gov/faqs/what-naming-convention-landsat-collections-level-1-scenes

[LS_578]

https://github.com/dgketchum/Landsat578#-1

[ExampleLandSat]

https://console.cloud.google.com/storage/browser/gcp-public-data-landsat/LC08/01/044/034/LC08_L1GT_044034_20130330_20170310_01_T2?_ga=2.210779154.665659046.1615242530-37570621.1615242530

[LandSatSuffixFormat]

https://prd-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/atoms/files/LSDS-750_Landsat8_Level-0-Reformatted_DataFormatControlBook-v15.pdf (page 26 / 99)

[LandsatProcLevels]

https://www.usgs.gov/core-science-systems/nli/landsat/landsat-levels-processing

[LandsatL2Names]

https://www.usgs.gov/faqs/what-naming-convention-landsat-collection-2-level-1-and-level-2-scenes?qt-news_science_products=0#qt-news_science_products

[LandSatARDDocs]

https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/LSDS-1873_US-Landsat%20C1-ARD-DFCB-v7.pdf

geowatch.gis.geotiff.walk_geotiff_products(dpath, with_product_dirs=True, with_loose_images=True, recursive=True)[source]¶

Walks a file path and returns directories and files that look like standalone geotiff products.

Parameters:: dpath (str) – directory to search
Yields:: str – paths of files or recognized geotiff product bundle directories

Example

>>> # xdoctest: +REQUIRES(env:SLOW_DOCTEST)
>>> # xdoctest: +REQUIRES(--network)
>>> # Test on real landsat data
>>> from geowatch.gis.geotiff import *  # NOQA
>>> import geowatch
>>> from os.path import dirname
>>> product = geowatch.demo.landsat_demodata.grab_landsat_product()
>>> dpath = dirname(dirname(ub.peek(product['bands'])))
>>> print(list(walk_geotiff_products(dpath)))
>>> print(list(walk_geotiff_products(dpath, with_product_dirs=False)))