gdal-2

GDAL

Note: This page is generated from a Jupyter Notebook. In Jupyter, to run a bash command, you need to prefix it with !. Hence you see that prefixing all gdal commands (For instance: !gdalinfo). However, when using GDAL as a CLI in your terminal, you do not have to prefix it with !.

In [6]:
cd /Users/abharathi/Documents/gis_data/gdal-tools/

/Users/abharathi/Documents/gis_data/gdal-tools

Verifying GDAL is installed correctly

In [8]:
!gdalinfo --version

GDAL 3.5.2, released 2022/09/02

gdalinfo - info about GIS datasets

Use to get basic info about any spatial dataset. Combine with -stats option to get stats about raster layers.

Getting info about DEM files using gdalinfo

Sample srtm folder has 4 DEM files of Mt. Everest.

In [15]:
ls -lh srtm/

total 202632
-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   430B Nov 23 13:38 N28E086.hgt.aux.xml
-rw-r--r--@ 1 abharathi  staff    25M Sep  7  2020 N28E087.hgt
In [16]:
!gdalinfo srtm/N27E086.hgt

Driver: SRTMHGT/SRTMHGT File Format
Files: srtm/N27E086.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,28.000138888888888)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (  85.9998611,  28.0001389) ( 85d59'59.50"E, 28d 0' 0.50"N)
Lower Left  (  85.9998611,  26.9998611) ( 85d59'59.50"E, 26d59'59.50"N)
Upper Right (  87.0001389,  28.0001389) ( 87d 0' 0.50"E, 28d 0' 0.50"N)
Lower Right (  87.0001389,  26.9998611) ( 87d 0' 0.50"E, 26d59'59.50"N)
Center      (  86.5000000,  27.5000000) ( 86d30' 0.00"E, 27d30' 0.00"N)
Band 1 Block=3601x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768
  Unit Type: m
In [17]:
!gdalinfo -stats srtm/N28E086.hgt

Driver: SRTMHGT/SRTMHGT File Format
Files: srtm/N28E086.hgt
       srtm/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
  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

gdalbuildvrt to build virtual raster

Combine all 4 DEM files into a single virtual raster. This .vrt file is a text file with pointers to actual data files. Gdal allows us to treat this as a physical file and run commands against it.

In [19]:
ls srtm/*.hgt

srtm/N27E086.hgt  srtm/N27E087.hgt  srtm/N28E086.hgt  srtm/N28E087.hgt
In [26]:
# gdalbuildvrt <output> <source_files>
!gdalbuildvrt srtm/merged.vrt srtm/*.hgt

0...10...20...30...40...50...60...70...80...90...100 - done.
In [27]:
ls srtm/

N27E086.hgt          N28E086.hgt          N28E087.hgt
N27E087.hgt          N28E086.hgt.aux.xml  merged.vrt
In [25]:
# preview the virtual raster file
cat srtm/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>

Finding highest point in Mr. Everest.

In [28]:
!gdalinfo -stats srtm/merged.vrt

Driver: VRT/Virtual Raster
Files: srtm/merged.vrt
       srtm/N27E086.hgt
       srtm/N27E087.hgt
       srtm/N28E086.hgt
       srtm/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

gdal_translate to convert file formats

Convert between hundreds of formats supported by gdal. Here we convert the virtual raster to a physical file on disk.

In [29]:
#gdal_translate -of <output_format> <input> <output>
!gdal_translate -of GTiff srtm/merged.vrt srtm/merged.tif

Input file size is 7201, 7201
0...10...20...30...40...50...60...70...80...90...100 - done.
In [30]:
ls -lh srtm/merged.tif

-rw-r--r--  1 abharathi  staff    99M Nov 23 14:17 srtm/merged.tif

Compressing output when using gdal_translate

By default, there is no compression. Gdal supports several compression algorithms and can be specified using the -co arg. Some options include

In [33]:
%%time
!gdal_translate -of GTiff srtm/merged.vrt srtm/merged.tif -co COMPRESS=DEFLATE

Input file size is 7201, 7201
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 31.2 ms, sys: 27.3 ms, total: 58.5 ms
Wall time: 3.57 s
In [34]:
ls -lh srtm/merged.tif

-rw-r--r--  1 abharathi  staff    72M Nov 23 14:35 srtm/merged.tif

Now the filesize has dropped from 99Mb to 72Mb. About 28% savings. Enhance compression even more by correlating values of neighboring cells (PREDICTOR) and by saving output in blocks (TILED).

In [35]:
%%time
!gdal_translate -of GTiff srtm/merged.vrt srtm/merged.tif \
    -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2

Input file size is 7201, 7201
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 60.3 ms, sys: 39.3 ms, total: 99.7 ms
Wall time: 10.4 s
In [36]:
ls -lh srtm/merged.tif

-rw-r--r--  1 abharathi  staff    39M Nov 23 14:38 srtm/merged.tif
In [ ]: