GDAL/OGR Quickstart
This Notebook is derived from the original GDAL-OGR quickstart adapted to run interactively in a IPython notebook and is devided in two parts GDAL (raster data) and OGR (vector data)
gdal_translate
gdal_warp
or gdal_merge.py
gdaltindex
Get to know GDAL
You will find the demo data at /usr/local/share/data
. We want to have a look at the naturalearth dataset data in this quickstart. We want to work with a copy of the data. So the first step is to copy the data to your home directory.
cd /home/user
cp -R /usr/local/share/data/natural_earth2/ ./gdal_natural_earth
You will then find a series of geotiff file:
ls /home/user/gdal_natural_earth/HYP_50M_SR_W.*
The web browser is not happy with tif files so to have a quick look at the image we convert it to a web friendly format (PNG) using the imagemagic utility convert
(this requires few seconds)
convert /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/HYP_50M_SR_W.png
display < /home/user/gdal_natural_earth/HYP_50M_SR_W.png
Get information about the raster data with gdalinfo
gdalinfo /home/user/gdal_natural_earth/HYP_50M_SR_W.tif
Note:
Simple Format Translation
First get to know your drivers. The --formats
commandline switch of gdalinfo can be used to see a list of available format drivers.
Each format reports if it is:
gdalinfo --formats
gdalinfo --format GTiff
Translation
Translations are accomplished with the gdal_translate command. The default output format is GeoTIFF. The -of
flag is used to select an output format and the -co
flag is used to specify a creation option:
gdal_translate -of JPEG -co QUALITY=10 /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/HYP_50M_SR_W.jpg
display < /home/user/gdal_natural_earth/HYP_50M_SR_W.jpg
The -ot switch can be used to alter the output data type.
gdal_translate -ot Int16 /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/HYP_50M_SR_W_Int16.tif
Use gdalinfo to verify data type.
gdalinfo /home/user/gdal_natural_earth/HYP_50M_SR_W.tif | grep Band
gdalinfo /home/user/gdal_natural_earth/HYP_50M_SR_W_Int16.tif | grep Band
Resizing
The -outsize switch can be used to set the size of the output file.
gdal_translate -outsize 50% 50% /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/HYP_50M_SR_W_small.tif
Use gdalinfo to verify the size.
gdalinfo /home/user/gdal_natural_earth/HYP_50M_SR_W.tif | grep Size
gdalinfo /home/user/gdal_natural_earth/HYP_50M_SR_W_small.tif | grep Size
Rescaling
The -scale switch can be used to rescale the data range of a given image. Explicit control of the input and output ranges is also available. The gdalinfo -mm
switch can be used to see pixel min/max values.
The output will display a computed Min/Max value for each band in the raster imput (3 band in our case)
gdalinfo -mm /home/user/gdal_natural_earth/HYP_50M_SR_W_small.tif | grep Min/Max
To rescale a specific band we can use the "-scale_bn" syntax where bn is a band number (e.g. "-scale_2" for the 2nd band of the output dataset), in the example below we will rescale the 3 bands of the HYP_50M_SR_W GeoTiff to be in the range 0-255 :
gdal_translate -scale_1 62.000 255.000 0 255.000 \
-scale_2 85.000 255.000 0 255.000 \
-scale_3 79.000 255.000 0 255.000 \
/home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/HYP_50M_SR_W_scaled.tif
Check the results with gdalinfo
:
gdalinfo -mm /home/user/gdal_natural_earth/HYP_50M_SR_W_scaled.tif | grep Min/Max
Let's convert the scaled output to JPG for a quick display, notice the color rearrangement.
gdal_translate -of JPEG -co QUALITY=40 /home/user/gdal_natural_earth/HYP_50M_SR_W_scaled.tif /home/user/gdal_natural_earth/HYP_50M_SR_W_scaled.jpg
display < /home/user/gdal_natural_earth/HYP_50M_SR_W_scaled.jpg
Splitting
Let’s split our image into two with -srcwin
which makes a copy based on pixel/line location (xoff yoff xsize ysize). You also could use -projwin
and define the corners in georeferenced coordinates (ulx uly lrx lry).
gdal_translate -srcwin 0 0 5400 5400 /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/west.tif
gdal_translate -srcwin 5400 0 5400 5400 /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/east.tif
gdal_translate -of JPEG -co QUALITY=40 /home/user/gdal_natural_earth/east.tif /home/user/gdal_natural_earth/east.jpg
gdal_translate -of JPEG -co QUALITY=40 /home/user/gdal_natural_earth/west.tif /home/user/gdal_natural_earth/west.jpg
Raster tileindex with gdaltindex
You can build a shapefile as a raster tileindex. For every image a polygon is generated with the bounds of the extent of the polygon and the path to the file.
gdaltindex /home/user/gdal_natural_earth/index_natural_earth.shp /home/user/gdal_natural_earth/*st.tif
The command above just created a new ESRI Shape File (default vector format), we will see how to work on vector files later on in the OGR section.
Reprojecting
For this process we assume that HYP_50M_SR_W.tif has been properly created with bounds. As we saw before with gdalinfo HYP_50M_SR_W has a proper coordinate system set.
If the tif file we want work on doesn't have proper projection information (wich is the case in most .tif files when associated with a world file .tfw) it is possible to assign a coordinate system to the image with gdal_translate
and the flag -a_srs
e.g :
gdal_translate -a_srs WGS84 HYP_50M_SR_W.tif HYP_50M_SR_W_4326.tif
Given a proper georeferenced raster file the gdalwarp
command can be used to assign a new Spatiale Reference System to it. Here we reproject the WGS84 geographic image to the Mercator projection:
gdalwarp -t_srs '+proj=merc +datum=WGS84' /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/mercator.tif
Creating output file that is 7292P x 9625L. Processing input file /home/user/gdal_natural_earth/HYP_50M_SR_W.tif. 0...10...20...30...40...50...60...70...80...90...100 - done.
convert /home/user/gdal_natural_earth/mercator.tif /home/user/gdal_natural_earth/mercator.png
0
display < /home/user/gdal_natural_earth/mercator.png
OGR
Like we did with raster using gdalinfo
, is possible to retrieve descriptive information from a vector datasource using the command line too ogrinfo
.
Let's run ogrinfo
on the previously generated shapefile index_natural_earth.shp
:
ogrinfo /home/user/gdal_natural_earth/index_natural_earth.shp index_natural_earth
INFO: Open of `/home/user/gdal_natural_earth/index_natural_earth.shp' using driver `ESRI Shapefile' successful. Layer name: index_natural_earth Geometry: Polygon Feature Count: 2 Extent: (-180.000000, -90.000000) - (180.000000, 90.000000) Layer SRS WKT: GEOGCS["GCS_WGS_1984", DATUM["WGS_1984", SPHEROID["WGS_84",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]] location: String (254.0) OGRFeature(index_natural_earth):0 location (String) = /home/user/gdal_natural_earth/east.tif POLYGON ((-0.00000000001796 90.0,179.999999999964047 90.0,179.999999999964047 -89.999999999982009,-0.00000000001796 -89.999999999982009,-0.00000000001796 90.0)) OGRFeature(index_natural_earth):1 location (String) = /home/user/gdal_natural_earth/west.tif POLYGON ((-180.0 90.0,-0.00000000001796 90.0,-0.00000000001796 -89.999999999982009,-180.0 -89.999999999982009,-180.0 90.0))
The shape file we just created can not be rendered directly in the notebook. Later in the python tutorial we'll learn how to use libraries like shapely to render directly vector data.
For now to display the results of our processing in the notebook we'll use the shp2img
command line tool (provided by mapserver ).
shp2img
require a mapserver mapfile as input to generate a rendered image.
Below we'll write a mapfile with 3 layers :
Note: the shapefile color has been classified based on the shapefile attribute location
with a translarency to show the 2 raster images underneat.
echo 'MAP
EXTENT -180 -89.999999999982 179.999999999964 90
IMAGETYPE "png"
NAME "simplepolygon"
SIZE 600 600
STATUS ON
UNITS DD
OUTPUTFORMAT
NAME "png"
MIMETYPE "image/png"
DRIVER "AGG/PNG"
EXTENSION "png"
IMAGEMODE RGB
TRANSPARENT TRUE
END # OUTPUTFORMAT
PROJECTION
"proj=longlat"
"datum=WGS84"
"no_defs"
END # PROJECTION
LAYER
DATA "/home/user/gdal_natural_earth/west.tif"
EXTENT -180 -89.999999999982 -1.79596892913025e-11 90
METADATA
"ows_title" "west"
END # METADATA
NAME "west"
PROJECTION
"proj=longlat"
"datum=WGS84"
"no_defs"
END # PROJECTION
STATUS ON
TILEITEM "location"
TYPE RASTER
UNITS METERS
END # LAYER
LAYER
DATA "/home/user/gdal_natural_earth/east.tif"
EXTENT -1.79596892913025e-11 -89.999999999982 179.999999999964 90
METADATA
"ows_title" "east"
END # METADATA
NAME "east"
PROJECTION
"proj=longlat"
"datum=WGS84"
"no_defs"
END # PROJECTION
STATUS ON
TILEITEM "location"
TYPE RASTER
UNITS METERS
END # LAYER
LAYER
DATA "/home/user/gdal_natural_earth/index_natural_earth.shp"
EXTENT -180 -89.999999999982 179.999999999964 90
NAME "index_natural_earth"
PROJECTION
"proj=longlat"
"datum=WGS84"
"no_defs"
END # PROJECTION
STATUS ON
TILEITEM "location"
TYPE POLYGON
UNITS METERS
CLASS
NAME "/home/user/gdal_natural_earth/east.tif"
EXPRESSION ("[location]" ="/home/user/gdal_natural_earth/east.tif")
STYLE
OPACITY 50
COLOR 218 57 57
END # STYLE
STYLE
OUTLINECOLOR 0 0 0
WIDTH 0.26
END # STYLE
END # CLASS
CLASS
NAME "/home/user/gdal_natural_earth/west.tif"
EXPRESSION ("[location]" ="/home/user/gdal_natural_earth/west.tif")
STYLE
OPACITY 50
COLOR 18 211 14
END # STYLE
STYLE
OUTLINECOLOR 0 0 0
WIDTH 0.26
END # STYLE
END # CLASS
END # LAYER
END # MAP' > index_natural_earth.map
shp2img -m index_natural_earth.map -i PNG -o index_natural_earth.png
display < index_natural_earth.png
#gdalwarp -t_srs EPSG:4326 -te -180.0000000 -90.0000000 180.0000000 90.0000000 /home/user/gdal_natural_earth/HYP_50M_SR_W.tif /home/user/gdal_natural_earth/HYP_50M_SR_W_geo.tif
#gdalinfo /home/user/gdal_natural_earth/HYP_50M_SR_W_geo.tif