Skip to content

Instantly share code, notes, and snippets.

@phargogh
Created October 27, 2022 01:48
Show Gist options
  • Save phargogh/38e5523386eb04aa50202ed547e13323 to your computer and use it in GitHub Desktop.
Save phargogh/38e5523386eb04aa50202ed547e13323 to your computer and use it in GitHub Desktop.
Creating a GTiff with embedded ISO metadata
import numpy
from osgeo import gdal
from osgeo import osr
# Taken from Examples section of https://wiki.gcube-system.org/gcube/ISO_19115:2003/19139
METADATA = """\
<fileIdentifier>
<gco:CharacterString>33462e78-e5ab-11c3-737d-b3a61366d028</gco:CharacterString>
</fileIdentifier>
<language>
<gco:CharacterString>eng</gco:CharacterString>
</language>
<hierarchyLevel>
<MD_ScopeCode codeList="http://metadata.dgiwg.org/codelistRegistry?MD_ScopeCode" codeListValue="dataset"/>
</hierarchyLevel>
<contact>
<CI_ResponsibleParty>
<individualName>
<gco:CharacterString>Mario Rossi</gco:CharacterString>
</individualName>
<organisationName>
<gco:CharacterString>ISTI-CNR</gco:CharacterString>
</organisationName>
<contactInfo>
<CI_Contact>
<address>
<CI_Address>
<electronicMailAddress>
<gco:CharacterString>rossi@isti.cnr.it</gco:CharacterString>
</electronicMailAddress>
</CI_Address>
</address>
</CI_Contact>
</contactInfo>
<role>
<CI_RoleCode codeList="http://www.isotc211.org/2005/resources/codeList.xml#CI_RoleCode" codeListValue="pointOfContact">pointOfContact</CI_RoleCode>
</role>
</CI_ResponsibleParty>
</contact>
<dateStamp>
<gco:Date>2013-03-09</gco:Date>
</dateStamp>
<metadataStandardName>
<gco:CharacterString>ISO19119</gco:CharacterString>
</metadataStandardName>
<metadataStandardVersion>
<gco:CharacterString>2005/PDAM 1</gco:CharacterString>
</metadataStandardVersion>"""
origin = (444720, 3751320)
pixel_size = (30, -30)
epsg = 3116
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
projection_wkt = srs.ExportToWkt()
target_path = 'raster.tif'
base_array = numpy.ones((4, 4))
raster_driver = gdal.GetDriverByName('GTiff')
ny, nx = base_array.shape
new_raster = raster_driver.Create(
target_path, nx, ny, 1, gdal.GDT_Float32)
if projection_wkt is not None:
new_raster.SetProjection(projection_wkt)
new_raster.SetGeoTransform(
[origin[0], pixel_size[0], 0.0, origin[1], 0.0, pixel_size[1]])
# NOTE: this GEO_METADATA could really be whatever metadata tag we want.
# For example, "FOO_GEO_METADATA" would be written as expected, just undera nonstandard tag.
new_raster.SetMetadataItem('GEO_METADATA', METADATA)
new_band = new_raster.GetRasterBand(1)
new_band.SetNoDataValue(-1)
new_band.WriteArray(base_array)
new_band = None
pew_raster = None
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment