26. Geospatial Data Processing with GDAL and OGR#
26.1. Introduction#
26.1.1. When to Use GDAL Directly#
26.2. Learning Objectives#
26.3. Installation and Setup#
conda install -c conda-forge gdal pygis
26.4. Sample Datasets#
wget https://github.com/opengeos/datasets/releases/download/gdal/gdal_sample_data.zip
unzip gdal_sample_data.zip
26.5. Understanding Your Data#
26.5.1. Examining Raster Data#
gdalinfo dem.tif
26.5.2. Examining Vector Data#
ogrinfo buildings.geojson
ogrinfo buildings.geojson -al -so
26.6. Coordinate Transformation#
26.6.1. Reprojecting Raster Data#
gdalwarp -t_srs EPSG:4326 dem.tif dem_4326.tif
26.6.2. Reprojecting Vector Data#
ogr2ogr -t_srs EPSG:3857 buildings_3857.gpkg buildings.geojson
26.7. Format Conversion#
26.7.1. Converting Raster Formats#
gdal_translate dem.tif dem.img
26.7.2. Converting Vector Formats#
ogr2ogr buildings.gpkg buildings.geojson
ogr2ogr las_vegas_bbox.fgb las_vegas_bbox.geojson
26.8. Clipping and Masking#
26.8.1. Clipping Raster with Vector Boundaries#
gdalwarp -cutline las_vegas_bbox.geojson -crop_to_cutline landsat.tif las_vegas.tif
26.8.2. Clipping Vector Data#
ogr2ogr -clipsrc las_vegas_bbox.geojson las_vegas_roads_clipped.geojson las_vegas_roads.geojson
ogrinfo las_vegas_bbox_4326.geojson -al -so
ogr2ogr -clipsrc -115.387634 35.943333 -114.883495 36.359161 las_vegas_roads_clipped.geojson las_vegas_roads.geojson
26.9. Raster Analysis and Calculations#
26.9.1. Working with Individual Bands#
gdal_translate -b 5 landsat.tif nir.tif
gdal_translate -b 4 landsat.tif red.tif
gdal_translate -b 3 landsat.tif green.tif
26.9.2. Performing Band Mathematics#
gdal_calc.py \
-A nir.tif \
-B red.tif \
--outfile=ndvi.tif \
--calc="(A.astype(float) - B) / (A + B + 1e-6)" \
--type=Float32 \
--NoDataValue=-9999
gdal_calc.py \
-A ndvi.tif \
--outfile=ndvi_clipped.tif \
--calc="(A < -1)*-1 + (A > 1)*1 + ((A >= -1) * (A <= 1))*A" \
--type=Float32 \
--NoDataValue=-9999 \
--overwrite
gdal_calc.py -A ndvi_clipped.tif --outfile=vegetation.tif --calc="A>0.3"
gdal_translate -a_nodata 0 vegetation.tif vegetation_bin.tif
gdal_edit.py -a_nodata 0 vegetation.tif
26.10. Converting Between Raster and Vector#
26.10.1. Vectorization#
gdal_polygonize.py building_masks.tif building_masks.gpkg
26.10.2. Rasterization#
gdal_rasterize -a uid -tr 0.6 0.6 -l buildings buildings_3857.gpkg buildings.tif
gdal_rasterize -burn 1 -tr 0.6 0.6 -l buildings buildings_3857.gpkg buildings2.tif
26.11. Geometry Processing#
26.11.1. Simplifying Complex Geometries#
ogr2ogr -f GPKG -t_srs EPSG:26911 -simplify 1 simplified.gpkg building_masks.gpkg
26.11.2. Dissolving Features by Attributes#
ogrinfo -al -so us_states.geojson
ogr2ogr -dialect sqlite -sql "SELECT ST_Union(geometry), country FROM us_states GROUP BY country" us_states_dissolved.gpkg us_states.geojson
26.11.3. Exploding Multi-part Geometries#
ogr2ogr -explodecollections hawaii_parts.geojson hawaii.geojson
26.12. Managing Fields and Layers#
26.12.1. Selecting Specific Fields#
ogr2ogr -select id,name us_states_select.geojson us_states.geojson
26.12.2. Renaming Layers#
ogr2ogr -nln states us_states_rename.gpkg us_states_select.geojson
26.12.3. Adding New Fields#
ogrinfo us_states_rename.gpkg -sql "ALTER TABLE states ADD COLUMN area DOUBLE"
26.13. Tiling and Data Management#
26.13.1. Creating Raster Tiles#
mkdir -p tiles
gdal_retile.py -targetDir tiles -ps 512 512 -co "TILED=YES" landsat.tif
26.13.2. Merging Raster Data#
gdal_merge.py -o landsat_merged.tif -of GTiff tiles/*.tif
gdalwarp -of GTiff tiles/*.tif landsat_merged2.tif
26.13.3. Working with Multiple Vector Files#
mkdir -p states
ogrinfo -ro -al -geom=NO us_states.geojson | grep "id (" | awk '{print $4}' | while read id; do
ogr2ogr -f GeoJSON "states/${id}.geojson" us_states.geojson -where "id = '${id}'"
done
ogr2ogr -f "ESRI Shapefile" us_states_merged.shp us_states.geojson
ogr2ogr -f "ESRI Shapefile" -update -append \
us_states_merged.shp hawaii.geojson -nln us_states_merged
26.14. Advanced Raster Processing#
26.14.1. Resampling and Resolution Management#
gdalwarp -tr 100 100 -r average dem.tif dem_100m.tif
26.14.2. Creating Band Composites#
gdal_merge.py -separate -o composite.tif nir.tif red.tif green.tif
gdalbuildvrt -separate stack.vrt nir.tif red.tif green.tif
gdal_translate stack.vrt composite.tif -co COMPRESS=DEFLATE
26.14.3. Handling Missing Data#
gdal_fillnodata.py -md 5 -of GTiff dem.tif filled_dem.tif
gdal_translate -a_nodata 0 dem.tif dem_nodata.tif
26.14.4. Cloud Optimized GeoTIFF#
gdal_translate dem.tif dem_cog.tif -of COG -co COMPRESS=DEFLATE
26.15. Terrain Analysis#
26.15.1. Computing Slope#
gdaldem slope dem.tif slope.tif \
-of GTiff \
-compute_edges
gdaldem slope dem.tif slope_percent.tif \
-of GTiff \
-compute_edges \
-p
26.15.2. Computing Aspect#
gdaldem aspect dem.tif aspect.tif \
-of GTiff \
-compute_edges
26.15.3. Creating Hillshade Visualizations#
gdaldem hillshade dem.tif hillshade.tif
gdaldem hillshade dem.tif multidirectional_hillshade.tif -multidirectional
26.15.4. Creating Color Relief Maps#
cat << EOF > colormap.txt
500, 51, 51, 153
1000, 3, 147, 249
1500, 37, 211, 109
2000, 181, 240, 138
2500, 218, 208, 133
3000, 146, 115, 94
4000, 183, 163, 159
5000, 255, 255, 255
EOF
gdaldem color-relief dem.tif colormap.txt color_relief.tif
26.15.5. Combining Hillshade and Color Relief#
gdal_calc.py \
-A color_relief.tif \
-B hillshade.tif \
--outfile=shaded_relief.tif \
--calc="((A.astype(float) * 0.6) + (B.astype(float) * 0.4))" \
--type=Byte
gdal_calc.py \
-A color_relief.tif \
-B hillshade.tif \
--outfile=color_hillshade.tif \
--calc="A * (B / 255.0)" \
--type=Byte
gdal_calc.py -A color_relief.tif --A_band=1 -B hillshade.tif --outfile=r.tif --calc="A*(B/255.0)" --type=Byte
gdal_calc.py -A color_relief.tif --A_band=2 -B hillshade.tif --outfile=g.tif --calc="A*(B/255.0)" --type=Byte
gdal_calc.py -A color_relief.tif --A_band=3 -B hillshade.tif --outfile=b.tif --calc="A*(B/255.0)" --type=Byte
gdal_merge.py -separate -o color_shaded_relief.tif r.tif g.tif b.tif
26.15.6. Generating Contour Lines#
gdal_contour -i 100 dem.tif contours.shp