This Quick Guide is based on the workshop in http://download.osgeo.org/gdal/workshop/foss4ge2015/workshop_gdal.pdf
This document is authored by Even Rouault and (C) Copyright Spatialys 2015. It is licenced under Creative Commons Attribution 3.0 (or any later version at the licensee choice) : https://creativecommons.org/licenses/by/3.0/.
Links to documentation in this tutorial used the versioned documentation for GDAL 1.11 (http://gdal.org/1.11/). For the documentation of the latest version (generally in development), remove the « 1.11/ »
Workshop data used and associated rights:
paris.tif : extract of OpenStreetMap. (C) OpenStreetMap contributors :
http://www.openstreetmap.org/copyrightworld.tif : from OSGeo-Live sample datane_10m_admin_0_countries.* and ne_10m_admin_1_states_provinces_shp.*: from
OSGeo -Live sample data (and originally from http://www.naturalearthdata.com,
public domain)wellington_west and wellington_east.png : derived from
https://data.linz.govt.nz/layer/1870-wellington-03m-rural-aerial-photos-2012-2013 ,
Licensed by Wellington Regional Council for reuse under the Creative Commons
Attribution 3.0 New Zealand licence (https://data.linz.govt.nz/license/attribution-3-0-
new-zealand/). For the purpose of this workshop, they have been post-processed to
decrease their resolution, reduce to 256 colors, add collars and convert to PNG.MK_30m.tif and ML_30m.tif : derived from https://data.linz.govt.nz/layer/1768-nz-
8m-digital-elevation-model-2012/ . License Creative Commons Attribution 3.0 New
Zealand (https://data.linz.govt.nz/license/attribution-3-0-new-zealand/). For the
purpose of this workshop, they have been post-processed to decrease their resolution.geomatrix.tif: from GDAL autotest suite. X/MIT Licensem2frac10bit.l1b : from GDAL extended data test suite, X/MIT License (http://download.osgeo.org/gdal/data/l1b/m2frac10bit.l1b)
In [1]:
; ls gdalworkshop/*
gdalinfo is the utility you will use all the time to discover metadata about a raster. This will also enable us to get a practical knowledge of most of the concepts of the GDAL data model.
Documentation of the gdalinfo utility : http://gdal.org/1.11/gdalinfo.html
In [2]:
; gdalinfo gdalworkshop/world.tif
Description of output :
In [3]:
using RasterIO
RasterIO.versioninfo("--version")
Out[3]:
In [4]:
raster = RasterIO.openraster("gdalworkshop/world.tif")
Out[4]:
In [5]:
filelist = RasterIO.getfilelist(raster.dataset)
Out[5]:
In [6]:
raster.width, raster.height, raster.nband
Out[6]:
In [7]:
RasterIO.wktprojection(raster)
Out[7]:
In [8]:
gt = RasterIO.geotransform(raster)
Out[8]:
In [9]:
print("origin: $((gt[1], gt[4]))")
In [10]:
print("pixel size: $((gt[2], gt[6]))")
In [11]:
print("upper left: $(RasterIO.applygeotransform(gt, 0.0, 0.0))")
In [12]:
print("upper right: $(RasterIO.applygeotransform(gt, Float64(raster.width), 0.0))")
In [13]:
print("lower left: $(RasterIO.applygeotransform(gt, 0.0, Float64(raster.height)))")
In [14]:
print("lower right: $(RasterIO.applygeotransform(gt, Float64(raster.width), Float64(raster.height)))")
In [15]:
print("center: $(RasterIO.applygeotransform(gt, raster.width/2, raster.height/2))")
In [16]:
RasterIO.getmetadata(raster.dataset)
Out[16]:
In [17]:
RasterIO.getmetadata(raster.dataset, "IMAGE_STRUCTURE")
Out[17]:
In [18]:
for i in 1:raster.nband
band = RasterIO.getrasterband(raster.dataset, i)
interp = RasterIO.getcolorinterp(band)
interp_name = RasterIO.colorinterpname(interp)
(w,h) = RasterIO.getblocksize(band)
println("Band $i, block size $((w,h)), color interp $interp_name")
for j in 1:RasterIO._overviewcount(band)
ovr_band = RasterIO.getoverview(band, j)
xsize = RasterIO._getrasterbandxsize(ovr_band)
ysize = RasterIO._getrasterbandysize(ovr_band)
println("Overview $j: $xsize x $ysize")
end
end
In [19]:
; gdalinfo gdalworkshop/geomatrix.tif
You can notice that contrary to the first example we no longer have a Origin and Pixel Size reported, but instead a GeoTransform.
The GeoTransform is the geo-transformation matrix, which is in mathematics described as an affine transformation from the coordinates in the pixel space (col,row) to the coordinates of the projected space (X,Y), with col and row starting from 0 for the upper-left pixel (in pixel space, not necessarily in projected space!). This is a series of 6 values gt[0], gt[1], … gt[5] such as :
X = gt[0] + col * gt[1] + row * gt[2]
Y = gt[3] + col * gt[4] + row * gt[5]
This is exactly the computation done by gdal.ApplyGeoTransform:
[X,Y]=gdal.ApplyGeoTransform(gt,row,col)
When the gt[2] and gt[4] terms are non zero, the image is no longer « north-up » with respect to the projected space, and you may have rotation and/or shearing.
To generate a more familiar image with those rotation/shearing being applied, you can use
gdalwarp.
The related mathematics are detailed at http://en.wikipedia.org/wiki/Transformation_matrix
In [20]:
; gdalinfo gdalworkshop/m2frac10bit.l1b | head -n 50
Lots of information here! The main novelty is the report of GCPs, which stand for Ground Control Points, also called Tiepoints or Control Points. Those are points of the raster for which the geolocation is known.
Let's look at the first one :
GCP[ 0]: Id=, Info=
(2023.5,221.5) -> (34.5374,59.0067,0)
Id, Info are text fields identifying the ground control points. They are rarely used in practice and often let to empty.
2023.5 and 221.5 are respectively the column and row in the raster for this GCP. As you can see, decimal values can be used. Here it means that the GCP is located at the center of the pixel (2023,221).
34.5374, 59.0067 and 0 are respectively the longitude/easting, latitude/northing and altitude of the GCP in the coordinate space. Here the reported GCP coordinate system is a geographic coordinate system, so they are longitude and latitude. The altitude is in most of the case 0.
You can also notice the lack of a geo-transform : the raster isn't projected yet. This can be done with gdalwarp
If the GCPs annoy you in the gdalinfo output, you can suppress them with the -nogcp option:
In [21]:
; gdalinfo -nogcp gdalworkshop/m2frac10bit.l1b
In [22]:
; gdalinfo -stats -nogcp -nomd gdalworkshop/m2frac10bit.l1b
The -stats option asks gdalinfo to compute the minimum, maximum, mean and standard deviation for each band. The -nogcp and -nomd options are just to suppress GCP and metadata from the output.
Now run
In [23]:
; gdalinfo -nogcp -nomd gdalworkshop/m2frac10bit.l1b
compare the output. You may notice that a m2frac10bit.l1b.aux.xml file is now mentioned, but the statistics are still reported. This file has been generated when computing the statistics. Let's display it
In [24]:
; cat gdalworkshop/m2frac10bit.l1b.aux.xml
You can see this file collects the statistics (as well as the other metadata). In case you are wondering what PAM in the above stands for, it is Persistant Auxiliary Metadata. The PAM .aux.xml format is something GDAL specific (and of course indirectly recognized by all software based on GDAL).
Caution: when a .aux.xml file exists, even if you specify -stats, gdalinfo will not recompute the statistics from the raster file but will use the ones stored in the .aux.xml. So in case the raster file would have been updated with new values, you may need to manually delete the .aux.xml file.
This is the utility to use to do format conversions, subsetting, format optimization, adding georeferencing, ...
Documentation of the gdal_translate utility: http://gdal.org/1.11/gdal_translate.html
Try:
In [25]:
; gdal_translate gdalworkshop/wellington_west.png gdalworkshop/wellington_west.gif
You get a warning here mentionning that the conversion has been done to TIFF and not to GIF as it was perhaps intended. GDAL is generally not sensitive to extensions, so you have to be explicit about the format, otherwise the default output format, which is GeoTiff will be used.
To get the list of supported formats, do:
In [26]:
; gdal_translate --formats
The letters between parenthesits after the format short name indicate the capabilities :
ro means read-onlyrw means read and one-time creation onlyrw+ means read, create and updatev signs means that it support virtual files, we will see that later.So the GIF driver support read and one-time creation, which is enough for gdal_translate. So let's correct our last attempt to specify the output format with -of :
In [27]:
; gdal_translate -of GIF gdalworkshop/wellington_west.png gdalworkshop/wellington_west.gif
In [28]:
; gdalinfo -noct gdalworkshop/wellington_west.png
We can see the presence of georeferenced coordinates (Origin and Pixel Size), that comes from wellington_west.wld, a World File. This is a 6 line text file that gives the information needed to build the 6-value geo-transformation matrix. Note that the order of values and their semantics are not exactly the same (World File use center of pixel convention), but GDAL will make the needed conversions.
But here we lack the Coordinate System. We must use external knowedge here to fill the gap. This dataset being from New Zealand and the coordinates being obviously not longitudes/latitudes, we could assume that this dataset is in the New Zealand Transverse Mercator (NTZM) map projection (for example by searching New Zealand in http://epsg.io/?q=new+zealand), and this is indeed the case. This coordinate system is codified by the EPSG (European Petroleum Survey Group) as being « EPSG:2193 »
We can see its definition with:
In [29]:
; gdalsrsinfo EPSG:2193
Now let us create a GeoTIFF from that:
In [30]:
; gdal_translate -a_srs EPSG:2193 gdalworkshop/wellington_west.png gdalworkshop/wellington_west.tif
Let's look at the result:
In [31]:
; gdalinfo -noct gdalworkshop/wellington_west.tif
And let's check that the coordinate system we select is consistant by using the OGR geocoding API:
In [32]:
; ogrinfo :memory: -q -sql "SELECT ogr_geocode('Wellington, New Zealand')"
Let's suppose we are only interested in a 4 km x 4 km area centered around Wellington city center.
We will convert first the above longitude, latitude in NZGD2000 with :
$ gdaltransform -s_srs WGS84 -t_srs EPSG:2193
174.7772239 -41.2887639
The command expects tuple of longitude/easting latitude/northing to be entered on the console and validated by Enter. The application may be terminated with Ctrl+Z.
Output :
1748816.1501341 5427663.38112408 0
So if we substract/add 2000 from this coordinate, we can do:
In [33]:
; gdal_translate -projwin 1746816 5429663 1750816 5425663 gdalworkshop/wellington_west.tif gdalworkshop/wellington_city.tif
Conventions are « -projwin ulx uly lrx lry » where :
ulx: X of upper-left corneruly: Y of upper-left cornerlrx: X of lower-right cornerlry: Y of lower-right corner
In [34]:
; gdalinfo -mm gdalworkshop/MK_30m.tif
The -mm option computes only the min/max values. We can also see a NoData value being advertized. The NoData value is a special value which means that pixel that have that value have no valid information. The NoData value is ignored in min/max, statistics or histogram computations.
We might want to generate a more convenient visual dataset from the DEM with:
In [35]:
; gdal_translate -outsize 50% 50% -ot Byte -of PNG gdalworkshop/MK_30m.tif gdalworkshop/MK_vis.png
-outsize 50% 50% is just to get a half-size reduction. Note that gdal_translate uses nearest neighbour sampling, which is not always appropriate. gdalwarp can be used for other resampling methods.-ot Byte: asks for data type conversion
In [36]:
using Images
In [37]:
img = imread("gdalworkshop/MK_vis.png")
Out[37]:
You can notice large zones saturated to white. Why so ? Because -ot does only data type conversion and clip values that would become out-of-range. So any elevation above 255 got clipped to 255. We also want to rescale the range [0,928.151] to [0,255]. We can do this by adding -scale
In [38]:
; gdal_translate -outsize 50% 50% -ot Byte -scale -of PNG gdalworkshop/MK_30m.tif gdalworkshop/MK_vis.png
You can control more precisely the rescaling with:
-scale src_min src_max
or
-scale src_min src_max dst_min dst_max
When not specifying dst_min and dst_max, 0 and 255 are selected. When not specifying src_min and src_max, the minimum/maximum values of the dataset are used.
You can also experiment with non-linear scaling with the -exponent value. For example try with -exponent 0.9. This will give an image slightly brighter. For reference, the destination value (dst_val) is computed from the (src_val) according to the following formula :
dst_val = dst_min + (dst_max – dst_min) * ((src_val-src_min)/(src_max – src_min)) ^ exponent
In [ ]: