gdal_t2

GDAL day 2¶ ¶

GDAL - Geospatial Data Abstration Library

Set up¶ ¶

Use - opengeo micromamba env. Move to gdal-tools folder.

In [2]:

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

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

In [3]:

ls

1870_southern-india.jpg*      landsat8/
5d16a93f1cf0f6000579ad2c.tif  london_1m_dsm/
batch.py*                     naip/
batch_parallel.py*            precipitation.gpkg
earth_at_night.jpg*           prism/
earthquakes/                  spatial_query.gpkg
gdal_stuff.qgz                srtm/
geonames/                     worldcities.csv*

Talk to Cloud systems via streaming¶ ¶

In [1]:

!gdalinfo /vsigs/spatialthoughts-public-data/viirs_ntl_2021_global.tif --config GS_NO_SIGN_REQUEST YES

Driver: GTiff/GeoTIFF
Files: /vsigs/spatialthoughts-public-data/viirs_ntl_2021_global.tif
Size is 80152, 80196
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["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]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20038000.000000000000000,20049000.000000000000000)
Pixel Size = (500.000000000000000,-500.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
  LAYOUT=COG
  PREDICTOR=2
Corner Coordinates:
Upper Left  (-20038000.000,20049000.000) (179d59'44.10"E, 85d 3'36.09"N)
Lower Left  (-20038000.000,-20049000.000) (179d59'44.10"E, 85d 3'36.09"S)
Upper Right (20038000.000,20049000.000) (179d59'44.10"W, 85d 3'36.09"N)
Lower Right (20038000.000,-20049000.000) (179d59'44.10"W, 85d 3'36.09"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  NoData Value=-9999
  Overviews: 40076x40098, 20038x20049, 10019x10024, 5009x5012, 2504x2506, 1252x1253, 626x626, 313x313

Working with aerial imagery¶ ¶

In [34]:

cd ../naip

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

In [35]:

ls

aoi.cpg*        output_2_1.jp2* output_3_8.jp2* output_5_7.jp2* output_7_6.jp2*
aoi.dbf*        output_2_2.jp2* output_4_1.jp2* output_5_8.jp2* output_7_7.jp2*
aoi.prj*        output_2_3.jp2* output_4_2.jp2* output_6_1.jp2* output_7_8.jp2*
aoi.qpj*        output_2_4.jp2* output_4_3.jp2* output_6_2.jp2* output_8_1.jp2*
aoi.shp*        output_2_5.jp2* output_4_4.jp2* output_6_3.jp2* output_8_2.jp2*
aoi.shx*        output_2_6.jp2* output_4_5.jp2* output_6_4.jp2* output_8_3.jp2*
filelist.txt    output_2_7.jp2* output_4_6.jp2* output_6_5.jp2* output_8_4.jp2*
output_1_1.jp2* output_2_8.jp2* output_4_7.jp2* output_6_6.jp2* output_8_5.jp2*
output_1_2.jp2* output_3_1.jp2* output_4_8.jp2* output_6_7.jp2* output_8_6.jp2*
output_1_3.jp2* output_3_2.jp2* output_5_1.jp2* output_6_8.jp2* output_8_7.jp2*
output_1_4.jp2* output_3_3.jp2* output_5_2.jp2* output_7_1.jp2* output_8_8.jp2*
output_1_5.jp2* output_3_4.jp2* output_5_3.jp2* output_7_2.jp2*
output_1_6.jp2* output_3_5.jp2* output_5_4.jp2* output_7_3.jp2*
output_1_7.jp2* output_3_6.jp2* output_5_5.jp2* output_7_4.jp2*
output_1_8.jp2* output_3_7.jp2* output_5_6.jp2* output_7_5.jp2*

Mosaicing aerial images¶ ¶

In [6]:

ls *.jp2 > filelist.txt

Alpha band is an additional band that simply states which cells have value and which have no-Data. It is an image mask. Alpha is useful when you compress raster using PNG or JP2 formats and reduce quantization to 8bit Uint8. With such low dynamic range, it becomes tricky to use 0 or 255 as no-data.

In [10]:

!gdalinfo output_1_1.jp2

Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: output_1_1.jp2
Size is 5000, 5000
Coordinate System is:
PROJCRS["NAD83 / UTM zone 10N",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            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",0,
            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["North America - between 126°W and 120°W - onshore and offshore. Canada - British Columbia; Northwest Territories; Yukon. United States (USA) - California; Oregon; Washington."],
        BBOX[30.54,-126,81.8,-119.99]],
    ID["EPSG",26910]]
Data axis to CRS axis mapping: 1,2
Origin = (545028.000000000000000,4249325.999946000054479)
Pixel Size = (0.600000000000000,-0.600000000600000)
Image Structure Metadata:
  INTERLEAVE=PIXEL
  COMPRESSION_REVERSIBILITY=LOSSY
Corner Coordinates:
Upper Left  (  545028.000, 4249326.000) (122d29' 3.80"W, 38d23'27.65"N)
Lower Left  (  545028.000, 4246326.000) (122d29' 4.49"W, 38d21'50.32"N)
Upper Right (  548028.000, 4249326.000) (122d27' 0.14"W, 38d23'27.08"N)
Lower Right (  548028.000, 4246326.000) (122d27' 0.87"W, 38d21'49.76"N)
Center      (  546528.000, 4247826.000) (122d28' 2.32"W, 38d22'38.71"N)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Undefined
  Overviews: 2500x2500, 1250x1250, 625x625
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Undefined
  Overviews: 2500x2500, 1250x1250, 625x625
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Undefined
  Overviews: 2500x2500, 1250x1250, 625x625
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000

Build virtual tile that mosaics the drone images using `gdalbuildvrt`¶ ¶

In [36]:

%time
!gdalbuildvrt -input_file_list filelist.txt naip.vrt

CPU times: user 4 µs, sys: 1e+03 ns, total: 5 µs
Wall time: 12.9 µs
0...10...20...30...40...50...60...70...80...90...100 - done.

Create overview file for visualizing `gdal_translate`¶ ¶

In [38]:

%%time
!gdal_translate -of JPEG -outsize 2% 2% naip.vrt naip_preview.jpg

Input file size is 40000, 40000
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 38.4 ms, sys: 24.2 ms, total: 62.6 ms
Wall time: 3.22 s

In [39]:

%%time
!gdaltindex -write_absolute_path index.json --optfile filelist.txt

Creating new index file...
CPU times: user 20 ms, sys: 15.3 ms, total: 35.4 ms
Wall time: 1.17 s

Processing satellite image¶ ¶

In [45]:

cd ../landsat8/

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

In [46]:

ls -lh

total 1845384
-rwxr-xr-x@ 1 abharathi  staff   110M Oct 12  2019 RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B2.tif*
-rwxr-xr-x@ 1 abharathi  staff   115M Oct 12  2019 RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif*
-rwxr-xr-x@ 1 abharathi  staff   115M Oct 12  2019 RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif*
-rwxr-xr-x@ 1 abharathi  staff   131M Oct 12  2019 RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.tif*
-rwxr-xr-x@ 1 abharathi  staff   430M Oct 12  2019 RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B8.tif*
-rw-r--r--  1 abharathi  staff   2.5K Apr 26 06:47 rgb.vrt

Use `gdalbuildvrt` to composite bands into a virtual raster¶ ¶

In [30]:

!gdalbuildvrt -o rgb.vrt -separate \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B2.tif

0...10...20...30...40...50...60...70...80...90...100 - done.

In [47]:

!gdalbuildvrt -o allbands.vrt -separate \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B2.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.tif \
  RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B8.tif

0...10...20...30...40...50...60...70...80...90...100 - done.

Use `gdaltranslate` to convert virtual to persisted file¶ ¶

In [41]:

%%time
!gdal_translate rgp.vrt rgb.tif -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE

ERROR 4: rgp.vrt: No such file or directory
CPU times: user 7.96 ms, sys: 9.33 ms, total: 17.3 ms
Wall time: 631 ms

Build a visualization product that applies histogram stretch¶ ¶

In [50]:

%%time
!gdal_translate -scale 0 0.3 0 255 -ot Byte rgb.vrt rgb_stretch.tif

Input file size is 7701, 7861
Warning 1: for band 1, nodata value has been clamped to 0, the original value being out of range.
Warning 1: for band 2, nodata value has been clamped to 0, the original value being out of range.
Warning 1: for band 3, nodata value has been clamped to 0, the original value being out of range.
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 41.6 ms, sys: 21.8 ms, total: 63.4 ms
Wall time: 4.17 s

Raster algebra using `gdal_calc.py`¶ ¶

In [52]:

!which gdal_calc.py

/Users/abharathi/micromamba/envs/opengeo/bin/gdal_calc.py

In [66]:

%run -i /Users/abharathi/micromamba/envs/opengeo/bin/gdal_calc.py

usage: gdal_calc.py [--help] --calc [expression ...] [-a [filename ...]]
                    [--a_band [n ...]] [-b [filename ...]] [--b_band [n ...]]
                    [-c [filename ...]] [--c_band [n ...]] [-d [filename ...]]
                    [--d_band [n ...]] [-e [filename ...]] [--e_band [n ...]]
                    [-f [filename ...]] [--f_band [n ...]] [-g [filename ...]]
                    [--g_band [n ...]] [-h [filename ...]] [--h_band [n ...]]
                    [-i [filename ...]] [--i_band [n ...]] [-j [filename ...]]
                    [--j_band [n ...]] [-k [filename ...]] [--k_band [n ...]]
                    [-l [filename ...]] [--l_band [n ...]] [-m [filename ...]]
                    [--m_band [n ...]] [-n [filename ...]] [--n_band [n ...]]
                    [-o [filename ...]] [--o_band [n ...]] [-p [filename ...]]
                    [--p_band [n ...]] [-q [filename ...]] [--q_band [n ...]]
                    [-r [filename ...]] [--r_band [n ...]] [-s [filename ...]]
                    [--s_band [n ...]] [-t [filename ...]] [--t_band [n ...]]
                    [-u [filename ...]] [--u_band [n ...]] [-v [filename ...]]
                    [--v_band [n ...]] [-w [filename ...]] [--w_band [n ...]]
                    [-x [filename ...]] [--x_band [n ...]] [-y [filename ...]]
                    [--y_band [n ...]] [-z [filename ...]] [--z_band [n ...]]
                    [-A [filename ...]] [--A_band [n ...]] [-B [filename ...]]
                    [--B_band [n ...]] [-C [filename ...]] [--C_band [n ...]]
                    [-D [filename ...]] [--D_band [n ...]] [-E [filename ...]]
                    [--E_band [n ...]] [-F [filename ...]] [--F_band [n ...]]
                    [-G [filename ...]] [--G_band [n ...]] [-H [filename ...]]
                    [--H_band [n ...]] [-I [filename ...]] [--I_band [n ...]]
                    [-J [filename ...]] [--J_band [n ...]] [-K [filename ...]]
                    [--K_band [n ...]] [-L [filename ...]] [--L_band [n ...]]
                    [-M [filename ...]] [--M_band [n ...]] [-N [filename ...]]
                    [--N_band [n ...]] [-O [filename ...]] [--O_band [n ...]]
                    [-P [filename ...]] [--P_band [n ...]] [-Q [filename ...]]
                    [--Q_band [n ...]] [-R [filename ...]] [--R_band [n ...]]
                    [-S [filename ...]] [--S_band [n ...]] [-T [filename ...]]
                    [--T_band [n ...]] [-U [filename ...]] [--U_band [n ...]]
                    [-V [filename ...]] [--V_band [n ...]] [-W [filename ...]]
                    [--W_band [n ...]] [-X [filename ...]] [--X_band [n ...]]
                    [-Y [filename ...]] [--Y_band [n ...]] [-Z [filename ...]]
                    [--Z_band [n ...]] --outfile filename
                    [--NoDataValue value] [--hideNoData] [--type datatype]
                    [--format gdal_format] [--creation-option option]
                    [--allBands [a-z, A-Z]] [--overwrite] [--debug] [--quiet]
                    [--color-table COLOR_TABLE]
                    [--extent {ignore,fail,union,intersect} | --projwin ulx uly lrx lry]
                    [--projectionCheck]
gdal_calc.py: error: the following arguments are required: --calc, --outfile

An exception has occurred, use %tb to see the full traceback.

SystemExit: 2

Performing NDVI¶ ¶

In [ ]:

%%time
%run -i /Users/abharathi/micromamba/envs/opengeo/bin/gdal_calc.py \
  -A RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.TIF \
  -B RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.TIF \
  --outfile ndvi.tif --calc="(A-B)/(A+B)" --NoDataValue=-999

Exercise - creating FCC viz product¶ ¶

In [69]:

%%time
!gdalbuildvrt -o nrg.vrt -separate \
RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B5.tif \
RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B4.tif \
RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B3.tif

0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 5.41 ms, sys: 9.25 ms, total: 14.7 ms
Wall time: 424 ms

In [73]:

%%time
!gdal_translate -of JPEG nrg.vrt nrg.jpg -co COMPRESS=JPG -scale 0 0.3 0 255 -ot Byte -outsize 10% 10%

Input file size is 7701, 7861
Warning 1: for band 1, nodata value has been clamped to 0, the original value being out of range.
Warning 1: for band 2, nodata value has been clamped to 0, the original value being out of range.
Warning 1: for band 3, nodata value has been clamped to 0, the original value being out of range.
Warning 6: driver JPEG does not support creation option COMPRESS
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 9.38 ms, sys: 11.1 ms, total: 20.4 ms
Wall time: 589 ms

Pansharpen¶ ¶

In [76]:

!which gdal_pansharpen.py

/Users/abharathi/micromamba/envs/opengeo/bin/gdal_pansharpen.py

In [78]:

%%time
%run -i /Users/abharathi/micromamba/envs/opengeo/bin/gdal_pansharpen.py \
RT_LC08_L1TP_137042_20190920_20190926_01_T1_2019-09-20_B8.tif \
rgb.vrt pansharpened.vrt -r bilinear -co COMPRESS=DEFLATE -co PHOTOMETRIC=RGB

CPU times: user 6.39 ms, sys: 9.63 ms, total: 16 ms
Wall time: 68.6 ms

In [80]:

%%time
!gdal_translate -of JPEG pansharpened.vrt pansharpened.jpg -scale 0 0.3 0 255 -ot Byte -outsize 40% 40%

Input file size is 15401, 15721
Warning 1: for band 1, nodata value has been clamped to 0, the original value being out of range.
Warning 1: for band 2, nodata value has been clamped to 0, the original value being out of range.
Warning 1: for band 3, nodata value has been clamped to 0, the original value being out of range.
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 97.1 ms, sys: 40.8 ms, total: 138 ms
Wall time: 12.8 s

Exercise - Lidar of UK viewshed analysis¶ ¶

STEPS

mosaic using buildvrt
assign SRS of 4326 and persist raster
reproject using gdaltransform
gdal_viewshed to produce the viewshed raster

In [83]:

cd ../london_1m_dsm/

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

In [84]:

ls *.asc > filelist.txt

In [85]:

%%time
!gdalbuildvrt  -input_file_list filelist.txt londondsm.vrt

0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 25.1 ms, sys: 17.7 ms, total: 42.8 ms
Wall time: 1.55 s

In [86]:

# No CRS in the vrt
!gdalinfo londondsm.vrt

Driver: VRT/Virtual Raster
Files: londondsm.vrt
       tq3080_DSM_1M.asc
       tq3081_DSM_1M.asc
       tq3082_DSM_1M.asc
       tq3083_DSM_1M.asc
       tq3084_DSM_1M.asc
       tq3180_DSM_1M.asc
       tq3181_DSM_1M.asc
       tq3182_DSM_1M.asc
       tq3183_DSM_1M.asc
       tq3184_DSM_1M.asc
       tq3280_DSM_1M.asc
       tq3281_DSM_1M.asc
       tq3282_DSM_1M.asc
       tq3283_DSM_1M.asc
       tq3284_DSM_1M.asc
       tq3380_DSM_1M.asc
       tq3381_DSM_1M.asc
       tq3382_DSM_1M.asc
       tq3383_DSM_1M.asc
       tq3384_DSM_1M.asc
       tq3480_DSM_1M.asc
       tq3481_DSM_1M.asc
       tq3482_DSM_1M.asc
       tq3483_DSM_1M.asc
       tq3484_DSM_1M.asc
Size is 5000, 5000
Origin = (530000.000000000000000,185000.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left  (  530000.000,  185000.000) 
Lower Left  (  530000.000,  180000.000) 
Upper Right (  535000.000,  185000.000) 
Lower Right (  535000.000,  180000.000) 
Center      (  532500.000,  182500.000) 
Band 1 Block=128x128 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999

In [89]:

%%time
#Assign SRS
!gdal_translate londondsm.vrt londondsm.tif \
  -co COMPRESS=DEFLATE -co PREDICTOR=2 \
  -a_srs EPSG:27700

Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
CPU times: user 30.6 ms, sys: 18.2 ms, total: 48.8 ms
Wall time: 3.3 s

Reproject¶ ¶