GDAL - Inspecting rasters¶
GDAL - Geospatial Data Abstration Library
Set up¶
Use - opengeo
micromamba env.
Move to gdal-tools
folder.
cd /Users/abharathi/Documents/gis_data/gdal-tools/
/Users/abharathi/Documents/gis_data/gdal-tools
ls
1870_southern-india.jpg* geonames/ prism/ batch.py* landsat8/ spatial_query.gpkg batch_parallel.py* london_1m_dsm/ srtm/ earth_at_night.jpg* naip/ worldcities.csv* earthquakes/ precipitation.gpkg
Intro¶
GDAL and OGR started as two different programs, hence have a different usage pattern or design.
Why use GDAL
- faster
- smaller, fits headless, remote execution
- FOSS
- great for server-side execution
- great for compression
Speed
- GDAL is fast. If we use Python to kick off GDAL as a CLI - it is still fast (
os.system
) - If we use Python bindings, they are much slower than using raw GDAL
Tools with GDAL¶
We get 2 sets of tools - raster tools that come with GDAL and vector tools that come with OGR.
3 Main tools
gdalinfo
- quick metadatagdal_translate
- convert formats, apply compressiongdalwarp
- resampling and reprojecting
Other tools
gdaladdo
- pyramid tiles for images for quicker GUI loadinggdaltindex
- tile index will give overviews of tilesgdalbuiltvrt
- for creating virtual rastersgdaldem
,gdalcounter
- elevationgdal_rasterize
,gdal_poligonize
- for conversion to and from raster <--> vectorgdal_grid
- interpolationsnearblack
- work with edge artifacts in rasters
Python bindings
gdal_merge.py
- merge rastersgdal_calc.py
- raster algebragdal_pansharpen.py
- for optical imagesgdal_retile.py
- tiling / chipping for deep learning
OGR
ogrinfo
- lists source and metadataogr2ogr
- format conversionogrmerge.py
- merge vector layers
!which gdalinfo
/Users/abharathi/micromamba/envs/opengeo/bin/gdalinfo
!gdalinfo --version
GDAL 3.6.2, released 2023/01/02
ls /Users/abharathi/micromamba/envs/opengeo/bin/ogr*
/Users/abharathi/micromamba/envs/opengeo/bin/ogr2ogr* /Users/abharathi/micromamba/envs/opengeo/bin/ogr_layer_algebra.py* /Users/abharathi/micromamba/envs/opengeo/bin/ogrinfo* /Users/abharathi/micromamba/envs/opengeo/bin/ogrlineref* /Users/abharathi/micromamba/envs/opengeo/bin/ogrmerge.py* /Users/abharathi/micromamba/envs/opengeo/bin/ogrtindex*
Exercises¶
# formats supported by gdal
!gdalinfo --formats
Supported Formats: VRT -raster,multidimensional raster- (rw+v): Virtual Raster DERIVED -raster- (ro): Derived datasets using VRT pixel functions GTiff -raster- (rw+vs): GeoTIFF COG -raster- (wv): Cloud optimized GeoTIFF generator NITF -raster- (rw+vs): National Imagery Transmission Format RPFTOC -raster- (rovs): Raster Product Format TOC format ECRGTOC -raster- (rovs): ECRG TOC format HFA -raster- (rw+v): Erdas Imagine Images (.img) SAR_CEOS -raster- (rov): CEOS SAR Image CEOS -raster- (rov): CEOS Image JAXAPALSAR -raster- (rov): JAXA PALSAR Product Reader (Level 1.1/1.5) GFF -raster- (rov): Ground-based SAR Applications Testbed File Format (.gff) ELAS -raster- (rw+v): ELAS ESRIC -raster- (rov): Esri Compact Cache AIG -raster- (rov): Arc/Info Binary Grid AAIGrid -raster- (rwv): Arc/Info ASCII Grid GRASSASCIIGrid -raster- (rov): GRASS ASCII Grid ISG -raster- (rov): International Service for the Geoid SDTS -raster- (rov): SDTS Raster DTED -raster- (rwv): DTED Elevation Raster PNG -raster- (rwv): Portable Network Graphics JPEG -raster- (rwv): JPEG JFIF MEM -raster,multidimensional raster- (rw+): In Memory Raster JDEM -raster- (rov): Japanese DEM (.mem) GIF -raster- (rwv): Graphics Interchange Format (.gif) BIGGIF -raster- (rov): Graphics Interchange Format (.gif) ESAT -raster- (rov): Envisat Image Format FITS -raster,vector- (rw+): Flexible Image Transport System BSB -raster- (rov): Maptech BSB Nautical Charts XPM -raster- (rwv): X11 PixMap Format BMP -raster- (rw+v): MS Windows Device Independent Bitmap DIMAP -raster- (rovs): SPOT DIMAP AirSAR -raster- (rov): AirSAR Polarimetric Image RS2 -raster- (rovs): RadarSat 2 XML Product SAFE -raster- (rov): Sentinel-1 SAR SAFE Product PCIDSK -raster,vector- (rw+v): PCIDSK Database File PCRaster -raster- (rw+): PCRaster Raster File ILWIS -raster- (rw+v): ILWIS Raster Map SGI -raster- (rw+v): SGI Image File Format 1.0 SRTMHGT -raster- (rwv): SRTMHGT File Format Leveller -raster- (rw+v): Leveller heightfield Terragen -raster- (rw+v): Terragen heightfield netCDF -raster,multidimensional raster,vector- (rw+s): Network Common Data Format HDF4 -raster,multidimensional raster- (ros): Hierarchical Data Format Release 4 HDF4Image -raster- (rw+): HDF4 Dataset ISIS3 -raster- (rw+v): USGS Astrogeology ISIS cube (Version 3) ISIS2 -raster- (rw+v): USGS Astrogeology ISIS cube (Version 2) PDS -raster- (rov): NASA Planetary Data System PDS4 -raster,vector- (rw+vs): NASA Planetary Data System 4 VICAR -raster,vector- (rw+v): MIPL VICAR file TIL -raster- (rov): EarthWatch .TIL ERS -raster- (rw+v): ERMapper .ers Labelled JP2OpenJPEG -raster,vector- (rwv): JPEG-2000 driver based on OpenJPEG library L1B -raster- (rovs): NOAA Polar Orbiter Level 1b Data Set FIT -raster- (rwv): FIT Image GRIB -raster,multidimensional raster- (rwv): GRIdded Binary (.grb, .grb2) RMF -raster- (rw+v): Raster Matrix Format WCS -raster- (rovs): OGC Web Coverage Service WMS -raster- (rwvs): OGC Web Map Service MSGN -raster- (rov): EUMETSAT Archive native (.nat) RST -raster- (rw+v): Idrisi Raster A.1 GSAG -raster- (rwv): Golden Software ASCII Grid (.grd) GSBG -raster- (rw+v): Golden Software Binary Grid (.grd) GS7BG -raster- (rw+v): Golden Software 7 Binary Grid (.grd) COSAR -raster- (rov): COSAR Annotated Binary Matrix (TerraSAR-X) TSX -raster- (rov): TerraSAR-X Product COASP -raster- (ro): DRDC COASP SAR Processor Raster R -raster- (rwv): R Object Data Store MAP -raster- (rov): OziExplorer .MAP KMLSUPEROVERLAY -raster- (rwv): Kml Super Overlay WEBP -raster- (rwv): WEBP PDF -raster,vector- (rw+vs): Geospatial PDF Rasterlite -raster- (rwvs): Rasterlite MBTiles -raster,vector- (rw+v): MBTiles PLMOSAIC -raster- (ro): Planet Labs Mosaics API CALS -raster- (rwv): CALS (Type 1) WMTS -raster- (rwv): OGC Web Map Tile Service SENTINEL2 -raster- (rovs): Sentinel 2 MRF -raster- (rw+v): Meta Raster Format TileDB -raster- (rw+vs): TileDB PNM -raster- (rw+v): Portable Pixmap Format (netpbm) DOQ1 -raster- (rov): USGS DOQ (Old Style) DOQ2 -raster- (rov): USGS DOQ (New Style) PAux -raster- (rw+v): PCI .aux Labelled MFF -raster- (rw+v): Vexcel MFF Raster MFF2 -raster- (rw+): Vexcel MFF2 (HKV) Raster GSC -raster- (rov): GSC Geogrid FAST -raster- (rov): EOSAT FAST Format BT -raster- (rw+v): VTP .bt (Binary Terrain) 1.3 Format LAN -raster- (rw+v): Erdas .LAN/.GIS CPG -raster- (rov): Convair PolGASP NDF -raster- (rov): NLAPS Data Format EIR -raster- (rov): Erdas Imagine Raw DIPEx -raster- (rov): DIPEx LCP -raster- (rwv): FARSITE v.4 Landscape File (.lcp) GTX -raster- (rw+v): NOAA Vertical Datum .GTX LOSLAS -raster- (rov): NADCON .los/.las Datum Grid Shift NTv2 -raster- (rw+vs): NTv2 Datum Grid Shift CTable2 -raster- (rw+v): CTable2 Datum Grid Shift ACE2 -raster- (rov): ACE2 SNODAS -raster- (rov): Snow Data Assimilation System KRO -raster- (rw+v): KOLOR Raw ROI_PAC -raster- (rw+v): ROI_PAC raster RRASTER -raster- (rw+v): R Raster BYN -raster- (rw+v): Natural Resources Canada's Geoid ARG -raster- (rwv): Azavea Raster Grid format RIK -raster- (rov): Swedish Grid RIK (.rik) USGSDEM -raster- (rwv): USGS Optional ASCII DEM (and CDED) GXF -raster- (rov): GeoSoft Grid Exchange Format KEA -raster- (rw+v): KEA Image Format (.kea) BAG -raster,multidimensional raster,vector- (rw+v): Bathymetry Attributed Grid HDF5 -raster,multidimensional raster- (rovs): Hierarchical Data Format Release 5 HDF5Image -raster- (rov): HDF5 Dataset NWT_GRD -raster- (rw+v): Northwood Numeric Grid Format .grd/.tab NWT_GRC -raster- (rov): Northwood Classified Grid Format .grc/.tab ADRG -raster- (rw+vs): ARC Digitized Raster Graphics SRP -raster- (rovs): Standard Raster Product (ASRP/USRP) BLX -raster- (rwv): Magellan topo (.blx) PostGISRaster -raster- (rws): PostGIS Raster driver SAGA -raster- (rw+v): SAGA GIS Binary Grid (.sdat, .sg-grd-z) XYZ -raster- (rwv): ASCII Gridded XYZ HF2 -raster- (rwv): HF2/HFZ heightfield raster OZI -raster- (rov): OziExplorer Image File CTG -raster- (rov): USGS LULC Composite Theme Grid ZMap -raster- (rwv): ZMap Plus Grid NGSGEOID -raster- (rov): NOAA NGS Geoid Height Grids IRIS -raster- (rov): IRIS data (.PPI, .CAPPi etc) PRF -raster- (rov): Racurs PHOTOMOD PRF EEDAI -raster- (ros): Earth Engine Data API Image DAAS -raster- (ro): Airbus DS Intelligence Data As A Service driver SIGDEM -raster- (rwv): Scaled Integer Gridded DEM .sigdem TGA -raster- (rov): TGA/TARGA Image File Format OGCAPI -raster,vector- (rov): OGCAPI STACTA -raster- (rovs): Spatio-Temporal Asset Catalog Tiled Assets STACIT -raster- (rovs): Spatio-Temporal Asset Catalog Items GPKG -raster,vector- (rw+vs): GeoPackage CAD -raster,vector- (rovs): AutoCAD Driver PLSCENES -raster,vector- (ro): Planet Labs Scenes API NGW -raster,vector- (rw+s): NextGIS Web GenBin -raster- (rov): Generic Binary (.hdr Labelled) ENVI -raster- (rw+v): ENVI .hdr Labelled EHdr -raster- (rw+v): ESRI .hdr Labelled ISCE -raster- (rw+v): ISCE raster Zarr -raster,multidimensional raster- (rw+vs): Zarr HTTP -raster,vector- (ro): HTTP Fetching Wrapper
Basic raster processing¶
cd srtm/
/Users/abharathi/Documents/gis_data/gdal-tools/srtm
ls -lh
total 202624 -rw-r--r--@ 1 abharathi staff 25M Sep 7 2020 N27E086.hgt -rw-r--r--@ 1 abharathi staff 25M Sep 7 2020 N27E087.hgt -rw-r--r--@ 1 abharathi staff 25M Sep 7 2020 N28E086.hgt -rw-r--r--@ 1 abharathi staff 25M Sep 7 2020 N28E087.hgt
!gdalinfo N28E086.hgt
Driver: SRTMHGT/SRTMHGT File Format Files: N28E086.hgt Size is 3601, 3601 Coordinate System is: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (85.999861111111116,29.000138888888888) Pixel Size = (0.000277777777778,-0.000277777777778) Metadata: AREA_OR_POINT=Point Corner Coordinates: Upper Left ( 85.9998611, 29.0001389) ( 85d59'59.50"E, 29d 0' 0.50"N) Lower Left ( 85.9998611, 27.9998611) ( 85d59'59.50"E, 27d59'59.50"N) Upper Right ( 87.0001389, 29.0001389) ( 87d 0' 0.50"E, 29d 0' 0.50"N) Lower Right ( 87.0001389, 27.9998611) ( 87d 0' 0.50"E, 27d59'59.50"N) Center ( 86.5000000, 28.5000000) ( 86d30' 0.00"E, 28d30' 0.00"N) Band 1 Block=3601x1 Type=Int16, ColorInterp=Undefined NoData Value=-32768 Unit Type: m
# Can chain command options
!gdalinfo -hist N28E086.hgt
Driver: SRTMHGT/SRTMHGT File Format Files: N28E086.hgt N28E086.hgt.aux.xml Size is 3601, 3601 Coordinate System is: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (85.999861111111116,29.000138888888888) Pixel Size = (0.000277777777778,-0.000277777777778) Metadata: AREA_OR_POINT=Point Corner Coordinates: Upper Left ( 85.9998611, 29.0001389) ( 85d59'59.50"E, 29d 0' 0.50"N) Lower Left ( 85.9998611, 27.9998611) ( 85d59'59.50"E, 27d59'59.50"N) Upper Right ( 87.0001389, 29.0001389) ( 87d 0' 0.50"E, 29d 0' 0.50"N) Lower Right ( 87.0001389, 27.9998611) ( 87d 0' 0.50"E, 27d59'59.50"N) Center ( 86.5000000, 28.5000000) ( 86d30' 0.00"E, 28d30' 0.00"N) Band 1 Block=3601x1 Type=Int16, ColorInterp=Undefined Min=2535.000 Max=8275.000 Minimum=2535.000, Maximum=8275.000, Mean=5067.325, StdDev=545.184 256 buckets from 2523.75 to 8286.25: 7 11 6 4 11 46 36 48 45 42 94 143 242 298 395 482 496 466 429 435 506 491 638 659 734 758 843 861 920 1014 1152 1228 1155 1346 1392 1670 1641 1890 1869 1964 1975 2201 2333 2765 2763 2959 2951 3295 3200 3673 3604 3937 3987 4476 4344 4862 4650 5086 5042 5889 5871 6668 6716 7342 7326 7837 8373 9414 9173 9968 10277 12213 11758 13969 15494 17512 19362 19328 57636 128726 140793 208505 241340 179825 185147 168703 163827 153862 162822 161140 173363 170846 180251 175639 188740 179526 186968 184531 196035 193185 206067 199861 211957 200760 212829 206817 220143 209265 222936 216445 229534 229816 238752 226566 234329 213829 218165 207523 217529 210291 219284 203776 206762 188260 198659 187907 184265 175861 162441 162943 150038 146503 129777 126385 111852 106263 94471 92537 81691 80049 71584 69721 62924 62460 56267 56311 51676 51844 46486 46417 42682 42311 39751 40405 37214 37815 35293 35475 33431 34307 32683 33092 30774 30803 28106 28282 25729 26323 24228 23774 22093 21084 18592 18101 16228 16004 14372 13961 12660 11036 10310 8965 8688 7677 7532 7126 6839 6015 5821 5071 4685 4092 3778 3257 3316 3134 2933 2527 2506 2178 2151 1975 1926 1736 1651 1552 1409 1222 1216 1070 1095 903 854 841 796 763 802 752 729 659 663 649 703 653 633 594 563 516 576 567 461 433 419 405 314 331 288 249 182 201 179 145 119 119 107 144 125 100 79 92 49 23 21 15 19 6 NoData Value=-32768 Unit Type: m Metadata: STATISTICS_MAXIMUM=8275 STATISTICS_MEAN=5067.3254290577 STATISTICS_MINIMUM=2535 STATISTICS_STDDEV=545.18375677872 STATISTICS_VALID_PERCENT=100
Making virtual rasters¶
Helps save space. A .vrt
file is a pointer to the source files and acts like a container.
ls *.hgt > srtm_filelist.txt
cat srtm_filelist.txt
N27E086.hgt N27E087.hgt N28E086.hgt N28E087.hgt
!gdalbuildvrt -input_file_list srtm_filelist.txt merged.vrt
0...10...20...30...40...50...60...70...80...90...100 - done.
cat merged.vrt
<VRTDataset rasterXSize="7201" rasterYSize="7201"> <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS> <GeoTransform> 8.5999861111111116e+01, 2.7777777777777778e-04, 0.0000000000000000e+00, 2.9000138888888888e+01, 0.0000000000000000e+00, -2.7777777777777778e-04</GeoTransform> <VRTRasterBand dataType="Int16" band="1"> <NoDataValue>-32768</NoDataValue> <ComplexSource> <SourceFilename relativeToVRT="1">N27E086.hgt</SourceFilename> <SourceBand>1</SourceBand> <SourceProperties RasterXSize="3601" RasterYSize="3601" DataType="Int16" BlockXSize="3601" BlockYSize="1" /> <SrcRect xOff="0" yOff="0" xSize="3601" ySize="3601" /> <DstRect xOff="0" yOff="3600" xSize="3601" ySize="3601" /> <NODATA>-32768</NODATA> </ComplexSource> <ComplexSource> <SourceFilename relativeToVRT="1">N27E087.hgt</SourceFilename> <SourceBand>1</SourceBand> <SourceProperties RasterXSize="3601" RasterYSize="3601" DataType="Int16" BlockXSize="3601" BlockYSize="1" /> <SrcRect xOff="0" yOff="0" xSize="3601" ySize="3601" /> <DstRect xOff="3600" yOff="3600" xSize="3601" ySize="3601" /> <NODATA>-32768</NODATA> </ComplexSource> <ComplexSource> <SourceFilename relativeToVRT="1">N28E086.hgt</SourceFilename> <SourceBand>1</SourceBand> <SourceProperties RasterXSize="3601" RasterYSize="3601" DataType="Int16" BlockXSize="3601" BlockYSize="1" /> <SrcRect xOff="0" yOff="0" xSize="3601" ySize="3601" /> <DstRect xOff="0" yOff="0" xSize="3601" ySize="3601" /> <NODATA>-32768</NODATA> </ComplexSource> <ComplexSource> <SourceFilename relativeToVRT="1">N28E087.hgt</SourceFilename> <SourceBand>1</SourceBand> <SourceProperties RasterXSize="3601" RasterYSize="3601" DataType="Int16" BlockXSize="3601" BlockYSize="1" /> <SrcRect xOff="0" yOff="0" xSize="3601" ySize="3601" /> <DstRect xOff="3600" yOff="0" xSize="3601" ySize="3601" /> <NODATA>-32768</NODATA> </ComplexSource> </VRTRasterBand> </VRTDataset>
!gdalinfo -stats merged.vrt
Driver: VRT/Virtual Raster Files: merged.vrt N27E086.hgt N27E087.hgt N28E086.hgt N28E087.hgt Size is 7201, 7201 Coordinate System is: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (85.999861111111116,29.000138888888888) Pixel Size = (0.000277777777778,-0.000277777777778) Corner Coordinates: Upper Left ( 85.9998611, 29.0001389) ( 85d59'59.50"E, 29d 0' 0.50"N) Lower Left ( 85.9998611, 26.9998611) ( 85d59'59.50"E, 26d59'59.50"N) Upper Right ( 88.0001389, 29.0001389) ( 88d 0' 0.50"E, 29d 0' 0.50"N) Lower Right ( 88.0001389, 26.9998611) ( 88d 0' 0.50"E, 26d59'59.50"N) Center ( 87.0000000, 28.0000000) ( 87d 0' 0.00"E, 28d 0' 0.00"N) Band 1 Block=128x128 Type=Int16, ColorInterp=Undefined Minimum=176.000, Maximum=8748.000, Mean=3840.966, StdDev=1649.758 NoData Value=-32768 Metadata: STATISTICS_MAXIMUM=8748 STATISTICS_MEAN=3840.9655621886 STATISTICS_MINIMUM=176 STATISTICS_STDDEV=1649.7576436136 STATISTICS_VALID_PERCENT=100
# Get maximum of the merged rasters
%time
!gdalinfo -stats -json merged.vrt | jq ".bands[0].maximum"
CPU times: user 2 µs, sys: 1e+03 ns, total: 3 µs Wall time: 9.06 µs 8748
Converting formats - gdal_translate
¶
!gdal_translate --formats
Supported Formats: VRT -raster,multidimensional raster- (rw+v): Virtual Raster DERIVED -raster- (ro): Derived datasets using VRT pixel functions GTiff -raster- (rw+vs): GeoTIFF COG -raster- (wv): Cloud optimized GeoTIFF generator NITF -raster- (rw+vs): National Imagery Transmission Format RPFTOC -raster- (rovs): Raster Product Format TOC format ECRGTOC -raster- (rovs): ECRG TOC format HFA -raster- (rw+v): Erdas Imagine Images (.img) SAR_CEOS -raster- (rov): CEOS SAR Image CEOS -raster- (rov): CEOS Image JAXAPALSAR -raster- (rov): JAXA PALSAR Product Reader (Level 1.1/1.5) GFF -raster- (rov): Ground-based SAR Applications Testbed File Format (.gff) ELAS -raster- (rw+v): ELAS ESRIC -raster- (rov): Esri Compact Cache AIG -raster- (rov): Arc/Info Binary Grid AAIGrid -raster- (rwv): Arc/Info ASCII Grid GRASSASCIIGrid -raster- (rov): GRASS ASCII Grid ISG -raster- (rov): International Service for the Geoid SDTS -raster- (rov): SDTS Raster DTED -raster- (rwv): DTED Elevation Raster PNG -raster- (rwv): Portable Network Graphics JPEG -raster- (rwv): JPEG JFIF MEM -raster,multidimensional raster- (rw+): In Memory Raster JDEM -raster- (rov): Japanese DEM (.mem) GIF -raster- (rwv): Graphics Interchange Format (.gif) BIGGIF -raster- (rov): Graphics Interchange Format (.gif) ESAT -raster- (rov): Envisat Image Format FITS -raster,vector- (rw+): Flexible Image Transport System BSB -raster- (rov): Maptech BSB Nautical Charts XPM -raster- (rwv): X11 PixMap Format BMP -raster- (rw+v): MS Windows Device Independent Bitmap DIMAP -raster- (rovs): SPOT DIMAP AirSAR -raster- (rov): AirSAR Polarimetric Image RS2 -raster- (rovs): RadarSat 2 XML Product SAFE -raster- (rov): Sentinel-1 SAR SAFE Product PCIDSK -raster,vector- (rw+v): PCIDSK Database File PCRaster -raster- (rw+): PCRaster Raster File ILWIS -raster- (rw+v): ILWIS Raster Map SGI -raster- (rw+v): SGI Image File Format 1.0 SRTMHGT -raster- (rwv): SRTMHGT File Format Leveller -raster- (rw+v): Leveller heightfield Terragen -raster- (rw+v): Terragen heightfield netCDF -raster,multidimensional raster,vector- (rw+s): Network Common Data Format HDF4 -raster,multidimensional raster- (ros): Hierarchical Data Format Release 4 HDF4Image -raster- (rw+): HDF4 Dataset ISIS3 -raster- (rw+v): USGS Astrogeology ISIS cube (Version 3) ISIS2 -raster- (rw+v): USGS Astrogeology ISIS cube (Version 2) PDS -raster- (rov): NASA Planetary Data System PDS4 -raster,vector- (rw+vs): NASA Planetary Data System 4 VICAR -raster,vector- (rw+v): MIPL VICAR file TIL -raster- (rov): EarthWatch .TIL ERS -raster- (rw+v): ERMapper .ers Labelled JP2OpenJPEG -raster,vector- (rwv): JPEG-2000 driver based on OpenJPEG library L1B -raster- (rovs): NOAA Polar Orbiter Level 1b Data Set FIT -raster- (rwv): FIT Image GRIB -raster,multidimensional raster- (rwv): GRIdded Binary (.grb, .grb2) RMF -raster- (rw+v): Raster Matrix Format WCS -raster- (rovs): OGC Web Coverage Service WMS -raster- (rwvs): OGC Web Map Service MSGN -raster- (rov): EUMETSAT Archive native (.nat) RST -raster- (rw+v): Idrisi Raster A.1 GSAG -raster- (rwv): Golden Software ASCII Grid (.grd) GSBG -raster- (rw+v): Golden Software Binary Grid (.grd) GS7BG -raster- (rw+v): Golden Software 7 Binary Grid (.grd) COSAR -raster- (rov): COSAR Annotated Binary Matrix (TerraSAR-X) TSX -raster- (rov): TerraSAR-X Product COASP -raster- (ro): DRDC COASP SAR Processor Raster R -raster- (rwv): R Object Data Store MAP -raster- (rov): OziExplorer .MAP KMLSUPEROVERLAY -raster- (rwv): Kml Super Overlay WEBP -raster- (rwv): WEBP PDF -raster,vector- (rw+vs): Geospatial PDF Rasterlite -raster- (rwvs): Rasterlite MBTiles -raster,vector- (rw+v): MBTiles PLMOSAIC -raster- (ro): Planet Labs Mosaics API CALS -raster- (rwv): CALS (Type 1) WMTS -raster- (rwv): OGC Web Map Tile Service SENTINEL2 -raster- (rovs): Sentinel 2 MRF -raster- (rw+v): Meta Raster Format TileDB -raster- (rw+vs): TileDB PNM -raster- (rw+v): Portable Pixmap Format (netpbm) DOQ1 -raster- (rov): USGS DOQ (Old Style) DOQ2 -raster- (rov): USGS DOQ (New Style) PAux -raster- (rw+v): PCI .aux Labelled MFF -raster- (rw+v): Vexcel MFF Raster MFF2 -raster- (rw+): Vexcel MFF2 (HKV) Raster GSC -raster- (rov): GSC Geogrid FAST -raster- (rov): EOSAT FAST Format BT -raster- (rw+v): VTP .bt (Binary Terrain) 1.3 Format LAN -raster- (rw+v): Erdas .LAN/.GIS CPG -raster- (rov): Convair PolGASP NDF -raster- (rov): NLAPS Data Format EIR -raster- (rov): Erdas Imagine Raw DIPEx -raster- (rov): DIPEx LCP -raster- (rwv): FARSITE v.4 Landscape File (.lcp) GTX -raster- (rw+v): NOAA Vertical Datum .GTX LOSLAS -raster- (rov): NADCON .los/.las Datum Grid Shift NTv2 -raster- (rw+vs): NTv2 Datum Grid Shift CTable2 -raster- (rw+v): CTable2 Datum Grid Shift ACE2 -raster- (rov): ACE2 SNODAS -raster- (rov): Snow Data Assimilation System KRO -raster- (rw+v): KOLOR Raw ROI_PAC -raster- (rw+v): ROI_PAC raster RRASTER -raster- (rw+v): R Raster BYN -raster- (rw+v): Natural Resources Canada's Geoid ARG -raster- (rwv): Azavea Raster Grid format RIK -raster- (rov): Swedish Grid RIK (.rik) USGSDEM -raster- (rwv): USGS Optional ASCII DEM (and CDED) GXF -raster- (rov): GeoSoft Grid Exchange Format KEA -raster- (rw+v): KEA Image Format (.kea) BAG -raster,multidimensional raster,vector- (rw+v): Bathymetry Attributed Grid HDF5 -raster,multidimensional raster- (rovs): Hierarchical Data Format Release 5 HDF5Image -raster- (rov): HDF5 Dataset NWT_GRD -raster- (rw+v): Northwood Numeric Grid Format .grd/.tab NWT_GRC -raster- (rov): Northwood Classified Grid Format .grc/.tab ADRG -raster- (rw+vs): ARC Digitized Raster Graphics SRP -raster- (rovs): Standard Raster Product (ASRP/USRP) BLX -raster- (rwv): Magellan topo (.blx) PostGISRaster -raster- (rws): PostGIS Raster driver SAGA -raster- (rw+v): SAGA GIS Binary Grid (.sdat, .sg-grd-z) XYZ -raster- (rwv): ASCII Gridded XYZ HF2 -raster- (rwv): HF2/HFZ heightfield raster OZI -raster- (rov): OziExplorer Image File CTG -raster- (rov): USGS LULC Composite Theme Grid ZMap -raster- (rwv): ZMap Plus Grid NGSGEOID -raster- (rov): NOAA NGS Geoid Height Grids IRIS -raster- (rov): IRIS data (.PPI, .CAPPi etc) PRF -raster- (rov): Racurs PHOTOMOD PRF EEDAI -raster- (ros): Earth Engine Data API Image DAAS -raster- (ro): Airbus DS Intelligence Data As A Service driver SIGDEM -raster- (rwv): Scaled Integer Gridded DEM .sigdem TGA -raster- (rov): TGA/TARGA Image File Format OGCAPI -raster,vector- (rov): OGCAPI STACTA -raster- (rovs): Spatio-Temporal Asset Catalog Tiled Assets STACIT -raster- (rovs): Spatio-Temporal Asset Catalog Items GPKG -raster,vector- (rw+vs): GeoPackage CAD -raster,vector- (rovs): AutoCAD Driver PLSCENES -raster,vector- (ro): Planet Labs Scenes API NGW -raster,vector- (rw+s): NextGIS Web GenBin -raster- (rov): Generic Binary (.hdr Labelled) ENVI -raster- (rw+v): ENVI .hdr Labelled EHdr -raster- (rw+v): ESRI .hdr Labelled ISCE -raster- (rw+v): ISCE raster Zarr -raster,multidimensional raster- (rw+vs): Zarr HTTP -raster,vector- (ro): HTTP Fetching Wrapper
%time
# gdal will guess the output format
!gdal_translate merged.vrt merged.tif
CPU times: user 3 µs, sys: 1e+03 ns, total: 4 µs Wall time: 7.87 µs Input file size is 7201, 7201 0...10...20...30...40...50...60...70...80...90...100 - done.
ls -lh merged*
-rw-r--r-- 1 abharathi staff 99M Apr 25 07:15 merged.tif -rw-r--r-- 1 abharathi staff 2.6K Apr 25 06:58 merged.vrt
%time
# gdal will guess the output format
!gdal_translate merged.vrt merged.zarr
CPU times: user 2 µs, sys: 1e+03 ns, total: 3 µs Wall time: 11 µs ERROR 1: Cannot guess driver for merged.zarr Output driver not found.
!gdalinfo merged.tif
Driver: GTiff/GeoTIFF Files: merged.tif Size is 7201, 7201 Coordinate System is: GEOGCRS["WGS 84", ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)"], MEMBER["World Geodetic System 1984 (G730)"], MEMBER["World Geodetic System 1984 (G873)"], MEMBER["World Geodetic System 1984 (G1150)"], MEMBER["World Geodetic System 1984 (G1674)"], MEMBER["World Geodetic System 1984 (G1762)"], MEMBER["World Geodetic System 1984 (G2139)"], ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], USAGE[ SCOPE["Horizontal component of 3D system."], AREA["World."], BBOX[-90,-180,90,180]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (85.999861111111116,29.000138888888888) Pixel Size = (0.000277777777778,-0.000277777777778) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 85.9998611, 29.0001389) ( 85d59'59.50"E, 29d 0' 0.50"N) Lower Left ( 85.9998611, 26.9998611) ( 85d59'59.50"E, 26d59'59.50"N) Upper Right ( 88.0001389, 29.0001389) ( 88d 0' 0.50"E, 29d 0' 0.50"N) Lower Right ( 88.0001389, 26.9998611) ( 88d 0' 0.50"E, 26d59'59.50"N) Center ( 87.0000000, 28.0000000) ( 87d 0' 0.00"E, 28d 0' 0.00"N) Band 1 Block=7201x1 Type=Int16, ColorInterp=Gray Min=176.000 Max=8748.000 Minimum=176.000, Maximum=8748.000, Mean=3840.966, StdDev=1649.758 NoData Value=-32768 Metadata: STATISTICS_MAXIMUM=8748 STATISTICS_MEAN=3840.9655621886 STATISTICS_MINIMUM=176 STATISTICS_STDDEV=1649.7576436136 STATISTICS_VALID_PERCENT=100
Convert format and compress¶
Use creation options (-co
) to specify compression when creating the resulting raser. Doc: https://gdal.org/drivers/raster/gtiff.html#creation-options
TILDED=YES
- stores tiled TIFF where neighborhood effect is used for compression
PREDICTOR=YES
- stores diff in values
%time
!gdal_translate merged.vrt merged.tif -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2
CPU times: user 3 µs, sys: 1e+03 ns, total: 4 µs Wall time: 9.78 µs Input file size is 7201, 7201 0...10...20...30...40...50...60...70...80...90...100 - done.
ls -lh *.tif
-rw-r--r-- 1 abharathi staff 39M Apr 25 07:31 merged.tif
!gdalinfo merged.tif
Driver: GTiff/GeoTIFF Files: merged.tif Size is 7201, 7201 Coordinate System is: GEOGCRS["WGS 84", ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)"], MEMBER["World Geodetic System 1984 (G730)"], MEMBER["World Geodetic System 1984 (G873)"], MEMBER["World Geodetic System 1984 (G1150)"], MEMBER["World Geodetic System 1984 (G1674)"], MEMBER["World Geodetic System 1984 (G1762)"], MEMBER["World Geodetic System 1984 (G2139)"], ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], USAGE[ SCOPE["Horizontal component of 3D system."], AREA["World."], BBOX[-90,-180,90,180]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (85.999861111111116,29.000138888888888) Pixel Size = (0.000277777777778,-0.000277777777778) Metadata: AREA_OR_POINT=Area Image Structure Metadata: COMPRESSION=DEFLATE INTERLEAVE=BAND PREDICTOR=2 Corner Coordinates: Upper Left ( 85.9998611, 29.0001389) ( 85d59'59.50"E, 29d 0' 0.50"N) Lower Left ( 85.9998611, 26.9998611) ( 85d59'59.50"E, 26d59'59.50"N) Upper Right ( 88.0001389, 29.0001389) ( 88d 0' 0.50"E, 29d 0' 0.50"N) Lower Right ( 88.0001389, 26.9998611) ( 88d 0' 0.50"E, 26d59'59.50"N) Center ( 87.0000000, 28.0000000) ( 87d 0' 0.00"E, 28d 0' 0.00"N) Band 1 Block=256x256 Type=Int16, ColorInterp=Gray Min=176.000 Max=8748.000 Minimum=176.000, Maximum=8748.000, Mean=3840.966, StdDev=1649.758 NoData Value=-32768 Metadata: STATISTICS_MAXIMUM=8748 STATISTICS_MEAN=3840.9655621886 STATISTICS_MINIMUM=176 STATISTICS_STDDEV=1649.7576436136 STATISTICS_VALID_PERCENT=100
Assign a different no data value¶
%time
!gdal_translate merged.vrt merged.tif -co COMPRESS=DEFLATE \
-co TILED=YES -co PREDICTOR=2 -a_nodata -9999
CPU times: user 4 µs, sys: 0 ns, total: 4 µs Wall time: 10 µs Input file size is 7201, 7201 0...10...20...30...40...50...60...70...80...90...100 - done.
!gdalinfo -json merged.tif | jq ".bands[0].noDataValue"
-9999
Writing GOG¶
Need to specify output format using -of COG
as COG also has same *.tif
extension
%time
!gdal_translate -of COG merged.vrt merged_cog.tif \
-co COMPRESS=DEFLATE -co PREDICTOR=2 -co NUM_THREADS=ALL_CPUS \
-a_nodata -9999
CPU times: user 4 µs, sys: 0 ns, total: 4 µs Wall time: 10 µs Input file size is 7201, 7201 0...10...20...30...40...50...60...70...80...90...100 - done.
!ls -lh *.tif
-rw-r--r-- 1 abharathi staff 39M Apr 25 07:40 merged.tif -rw-r--r-- 1 abharathi staff 54M Apr 25 07:50 merged_cog.tif
!gdalinfo merged_cog.tif -json | jq ".metadata.IMAGE_STRUCTURE"
{ "COMPRESSION": "DEFLATE", "INTERLEAVE": "BAND", "LAYOUT": "COG", "PREDICTOR": "2" }
Using GDALINFO on COG files stored on the cloud¶
I am using Open aerial map site for this: https://map.openaerialmap.org/#/-117.81875610351562,33.76544869849223,9/latest/643b06ca8cae390005a1456f?_k=9krcby
!gdalinfo https://oin-hotosm.s3.amazonaws.com/643b06548cae390005a1456c/0/643b06548cae390005a1456d.tif
Driver: GTiff/GeoTIFF Files: none associated Size is 12099, 9659 Coordinate System is: PROJCRS["WGS 84 / UTM zone 36S", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["UTM zone 36S", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",33, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",0.9996, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",500000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",10000000, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["(E)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["(N)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["Engineering survey, topographic mapping."], AREA["Between 30°E and 36°E, southern hemisphere between 80°S and equator, onshore and offshore. Burundi. Eswatini (Swaziland). Kenya. Malawi. Mozambique. Rwanda. South Africa. Tanzania. Uganda. Zambia. Zimbabwe."], BBOX[-80,30,0,36]], ID["EPSG",32736]] Data axis to CRS axis mapping: 1,2 Origin = (415959.438299251720309,7169051.803544469177723) Pixel Size = (0.019999931184024,-0.019998767072739) Metadata: AREA_OR_POINT=Area Image Structure Metadata: COMPRESSION=YCbCr JPEG INTERLEAVE=PIXEL JPEGTABLESMODE=1 JPEG_QUALITY=75 SOURCE_COLOR_SPACE=YCbCr Corner Coordinates: Upper Left ( 415959.438, 7169051.804) ( 32d 9'47.25"E, 25d35'37.10"S) Lower Left ( 415959.438, 7168858.635) ( 32d 9'47.20"E, 25d35'43.38"S) Upper Right ( 416201.417, 7169051.804) ( 32d 9'55.92"E, 25d35'37.15"S) Lower Right ( 416201.417, 7168858.635) ( 32d 9'55.88"E, 25d35'43.42"S) Center ( 416080.428, 7168955.219) ( 32d 9'51.56"E, 25d35'40.26"S) Band 1 Block=512x512 Type=Byte, ColorInterp=Red Description = red Min=0.000 Max=255.000 Minimum=0.000, Maximum=255.000, Mean=108.198, StdDev=76.230 Overviews: 6050x4830, 3025x2415, 1513x1208, 757x604, 379x302 Mask Flags: PER_DATASET Overviews of mask band: 6050x4830, 3025x2415, 1513x1208, 757x604, 379x302 Metadata: STATISTICS_MAXIMUM=255 STATISTICS_MEAN=108.19774015389 STATISTICS_MINIMUM=0 STATISTICS_STDDEV=76.230400100428 Band 2 Block=512x512 Type=Byte, ColorInterp=Green Description = green Min=0.000 Max=255.000 Minimum=0.000, Maximum=255.000, Mean=99.308, StdDev=66.015 Overviews: 6050x4830, 3025x2415, 1513x1208, 757x604, 379x302 Mask Flags: PER_DATASET Overviews of mask band: 6050x4830, 3025x2415, 1513x1208, 757x604, 379x302 Metadata: STATISTICS_MAXIMUM=255 STATISTICS_MEAN=99.30805629414 STATISTICS_MINIMUM=0 STATISTICS_STDDEV=66.014549775173 Band 3 Block=512x512 Type=Byte, ColorInterp=Blue Description = blue Min=0.000 Max=255.000 Minimum=0.000, Maximum=255.000, Mean=79.344, StdDev=57.837 Overviews: 6050x4830, 3025x2415, 1513x1208, 757x604, 379x302 Mask Flags: PER_DATASET Overviews of mask band: 6050x4830, 3025x2415, 1513x1208, 757x604, 379x302 Metadata: STATISTICS_MAXIMUM=255 STATISTICS_MEAN=79.343554107368 STATISTICS_MINIMUM=0 STATISTICS_STDDEV=57.837036128418
DEMs using gdaldem
¶
# We specify scale since the vert units are in meters whereas the horz units are in degrees
%time
!gdaldem hillshade merged.vrt hillshade.tif -s 111120
CPU times: user 4 µs, sys: 1 µs, total: 5 µs Wall time: 8.82 µs 0...10...20...30...40...50...60...70...80...90...100 - done.
ls -lh hill*
-rw-r--r-- 1 abharathi staff 49M Apr 25 08:00 hillshade.tif
DEM classified (color relief) file¶
# first create a file to store the classes
1000,101,146,82
1500,190,202,130
2000,241,225,145
2500,244,200,126
3000,197,147,117
4000,204,169,170
5000,251,238,253
6000,255,255,255
(6000, 255, 255, 255)
cat colormap.txt
1000,101,146,82 1500,190,202,130 2000,241,225,145 2500,244,200,126 3000,197,147,117 4000,204,169,170 5000,251,238,253 6000,255,255,255
%time
!gdaldem color-relief merged.vrt colormap.txt colorized_dem.tif
CPU times: user 2 µs, sys: 1e+03 ns, total: 3 µs Wall time: 11.2 µs 0...10...20...30...40...50...60...70...80...90...100 - done.
# Convert to PNG, but a downsampled file that is 10% x 10%
%time
!gdal_translate -of PNG colorized_dem.tif colorized.png -outsize 10% 10%
CPU times: user 5 µs, sys: 3 µs, total: 8 µs Wall time: 16.9 µs Input file size is 7201, 7201 0...10...20...30...40...50...60...70...80...90...100 - done.
ls -lh colorized*
-rw-r--r-- 1 abharathi staff 595K Apr 25 08:22 colorized.png -rw-r--r-- 1 abharathi staff 712B Apr 25 08:22 colorized.png.aux.xml -rw-r--r-- 1 abharathi staff 148M Apr 25 08:12 colorized_dem.tif
0.00009259259259 * 111320
10.3074074071188