Tutorial: Introductie GDAL

Site:	OpenCourseWare for GIS
Cursus:	Introductie GDAL & Python
Boek:	Tutorial: Introductie GDAL

Afgedrukt door:	Guest user
Datum:	zaterdag, 20 april 2024, 16:19

Inhoudsopgave

1. Introductie
2. Informatie opvragen uit GISdata
- 2.1. Informatie opvragen uit rasterdata
- 2.2. Informatie opvragen van vectordata
3. Herprojecteren van GIS data
- 3.1. Herprojecteer rasterdata
- 3.2. Herprojecteer vectordata
4. Converteren van GIS formaten
5. Batch conversie
6. Conclusie

1. Introductie

Wat is GDAL?

De Geospatial Data Abstraction Library (GDAL) is een open source computersoftwarebibliotheek voor het lezen en schrijven van raster- en vector gegevensformaten. Het is nuttig voor conversies van bestandsformaten en projecties. Er zijn ook tools voor GIS analyses en ruimtelijke queries. Zie voor meer info: https://gdal.org/.

Leerdoelen

Na deze cursus ben je in staat om:

Informatie op te vragen uit geografische data
Geografische data te herprojecteren
Rastereigenschappen te wijzigen
Rasterformaten te converteren
Vectorformaten te converteren
Ruimtelijke queries toe te passen op vectoren
Bestanden met door komma's gescheiden waarden (CSV bestanden) te converteren
Batchconversie uit te voeren

Software

Voor deze oefeningen moet GDAL worden geïnstalleerd, bij voorkeur met behulp van het OSGeo4W-distributiepakket. Als je QGIS hebt geïnstalleerd, is de OSGEO4W-distributie er al. In de video kun je zien hoe je de LTR-versie van QGIS installeert.

Data voor de oefeningen

De data voor de oefeningen kunnen worden gedownload van de hoofdpagina. Het bevat:

roads.shp: wegen van OpenStreetMap (http://openstreetmap.org)

srtm_37_02.tif: hoogtemodel van de Shuttle Radar Topography Mission (SRTM) (http://www2.jpl.nasa.gov/srtm/)

gem_2011_gn1.shp: grenzen van Nederlandse gemeented, open data van CBS en Kadaster (http://www.cbs.nl).

Locations.csv: tabel met lokaties van objecten.

landuse.zip: tijdserie van landgebruiksrasters in IDRISI formaat.

2. Informatie opvragen uit GISdata

In dit hoofdstuk leer je hoe je:

Informatie kan opvragen uit rasterdata
Informatie kan opvragen uit vectordata

We zullen de gdalinfo en ogrinfo commando's introduceren.

2.1. Informatie opvragen uit rasterdata

Een van de gemakkelijkste en meest bruikbare commando's in GDAL is gdalinfo. Wanneer een raster als argument wordt gegeven, wordt alle relevante informatie die bekend is over het bestand opgehaald en afgedrukt. Dit is vooral handig als het rasterbestand extra taggegevens bevat, zoals het geval is bij TIF-bestanden. Bij het werken met satellietbeelden is dit een uiterst handige manier om de locatie van de beelden in lat/lon coördinaten bij te houden, evenals de projectie informatie.

1. Open de opdrachtregelinterface waar je GDAL hebt geconfigureerd. Als je QGIS hebt geïnstalleerd kan je de OSGeo4W Shell gebruiken.

2. Ga met de opdrachtprompt naar de directory waar u de trainingsgegevens hebt opgeslagen, bijvoorbeeld Z:\gdal_exercises

3. Voer het volgende commando uit:

gdalinfo srtm_37_02.tif <ENTER>

4. Beantwoord de volgende vragen op basis van de informatie die op het scherm wordt getoond:

Wat is de grootte van het rasterbestand?
Wat is het coördinaatsysteem?
Wat is de EPSG code?
Wat zijn de minimum en maximum waarden?

Meer informatie over het gebruik van het gdalinfo commando kan je hier vinden: https://gdal.org/programs/gdalinfo.html

EPSG codes worden gebruikt om projecties te definiëren. Meer informatie hierover in de onderstaande video:

2.2. Informatie opvragen van vectordata

Om informatie van vectordata op te vragen gebruiken we het ogrinfo commando.

1. Voer het volgende commando uit:
ogrinfo -al roads.shp | more <ENTER>
ogrinfo -al gem_2011_gn1.shp | more <ENTER>

2. Beantwoord de volgende vragen:

Wat zijn de projecties van deze twee shapefiles?
Wat zijn de EPSG codes?

Zoals je kunt zien, zijn er verschillende EPSG-codes gerapporteerd. Dat komt omdat verschillende elementen van de projectie hun eigen EPSG-code hebben. Vuistregel is dat je de laatste EPSG-code kunt nemen.

Dus de projectie van roads.shp is EPSG: 4326, wat betekent dat het in het Geografisch Coordinaat Systeem (GCS) is met coordinaten in lengte- en breedtegraden.

Bij gem_2011_gn1.shp is de projectie EPSG: 28992. Dat is de Nederlandse projectie (Amersfoort RD/Nieuw).

Meer over het ogrinfo commando kan hier gevonden worden.

3. Herprojecteren van GIS data

In dit hoofdstuk leer je hoe je:

Rasterdata kan herprojecteren
Vectordata kan herprojecteren

We gebruiken hiervoor de gdalwarp en ogr2ogr commando's.

In deze oefening willen we een kaart maken van de gemeente Delft met de hoofdwegen en het reliëf. Omdat de datasets in verschillende formaten zijn, moeten we ze herprojecteren naar een gemeenschappelijk coördinatensysteem. Hier herprojecteren we alle datasets naar de Nederlandse Amersfoort/RD Nieuw projectie.

3.1. Herprojecteer rasterdata

GDAL heeft de mogelijkheid om een rastercoördinatensysteem te wijzigen met behulp van de volgende syntaxis:

gdalwarp -t_srs EPSG:... <input> <output>

Het -t_srs argument specificeert het doelcoördinatensysteem. Als het broncoördinatensysteem onbekend is, moet het worden opgegeven met het -s_srs argument. EPSG:... specificeert de EPSG-code van de projectie. <input> and <output> zijn respectievelijk de invoer- en uitvoergegevens.

We gaan nu een digitaal hoogtemodel van de Shuttle Radar Topography Mission (SRTM) herprojecteren. Je kan hoogtemodellen voor je eigen gebied downloaden van USGS Earth Explorer. Hier gebruiken we de data behorende bij deze cursus.

Om het hoogtemodel te herprojecteren van WGS-84 lat/lon naar Amersfoort/RD Nieuw gebruiken we dit commando:

gdalwarp -t_srs EPSG:XXXXX srtm_37_02.tif dem_rd.tif

1. Vervang XXXXX met de EPSG code voor Amersfoort/RD Nieuw (zie een van je eerdere antwoorden met ogrinfo).

2. Voer dit commando uit:

gdalwarp -t_srs EPSG:28992 srtm_37_02.tif dem_rd.tif <ENTER>

3. Visualiseer het resultaat in QGIS.

Meer informatie over het gdalwarp-commando is hier te vinden.

3.2. Herprojecteer vectordata

Voor vectordata wordt het ogr2ogr commando gebruikt. De syntax is:

ogr2ogr -t_srs EPSG:XXXXX <output> <input>

Merk op dat <output> en <input> in omgekeerde volgorde zijn dan bij gdalwarp! Dit verschil is algemeen bij gdal en ogr commando's.

Hier gaan we wegendata van OpenStreetMap converteren naar de Amersfoort/RD Nieuw projectie.

1. Voer het volgende commando uit:
ogr2ogr -t_srs EPSG:28992 roadsreprojected.shp roads.shp

Het commando wordt uitgevoerd (geef het wat tijd, het is een grote dataset), maar geeft een waarschuwing.

2. Wat betekent de waarschuwing?

Meer informatie over het ogr2ogr commando kan je hier vinden.

4. Converteren van GIS formaten

In dit hoofdstuk leer je hoe je:

rasterformaten kan converteren
vectorformaten kan converteren

We gebruiken hier de gdal_translate en ogr2ogr commando's.

4.1. Converteer rasterformaten

gdal_translate wordt gebruikt om rasterformaten te converteren. De algemene syntax is:

gdal_translate -of FORMAAT inputBestand outputBestand

Alle ondersteunde formaten kan je hier vinden.

We gaan nu het hoogtemodel van GeoTIFF formaat naar SAGA formaat converteren. SAGA is een GIS met z'n eigen, door GDAL ondersteunde, binaire formaat met de .sdat entensie.

1. Voer het volgende commando uit:
gdal_translate -of SAGA dem_rd.tif dem_rd.sdat

Sommige formaten hebben meer argumenten nodig. PCRaster bijvoorbeeld heeft een specificatie van het datatype nodig (boolean, nominal, scalar, etc.). Laten we dezelfde data naar PCRaster formaat converteren.

2. Voer het volgende commando uit:
gdal_translate -of PCRaster -ot Float32 -mo "VS_SCALAR" dem_rd.tif dem_rd.map

Om naar het PCRaster formaat the converteren moet je het datatype van het raster weten. In bovenstaand voorbeeld is het hoogtemodel een continu raster en heeft daarom het "scalar" datatype. Bij een discreet raster met klassen zou het datatype "nominal" zijn. Deze informatie is te vinden onder het rasterformaat op de GDAL website.

Data Type	-of	-mo
Boolean	Byte	"VS_BOOLEAN"
Nominal	Int32	"VS_NOMINAL"
Ordinal	Int32	"VS_ORDINAL"
Scalar	Float32 or Float64	"VS_SCALAR"
Direction	Float32 or Float64	"VS_DIRECTION"
LDD	Int32	"VS_LDD"

4.2. Converteer vectorformaten

Met het ogr2ogr commando kan je ook vectorformaten converteren.

Hier is een lijst met ondersteunde vectorformaten.

De algemene syntax is:

ogr2ogr -f <"Formaat"> <output> <input>

"Formaat" staat in de "short name" kolom in de vector drivers table en moet worden getypt tussen dubbele aanhalingstekens.

We gaan nu gem_2011_gn1.shp converteren naar een Google KML bestand dat geopend kan worden in Google Earth.

1. Voer het volgende commando uit:
ogr2ogr -f "KML" gem.kml gem_2011_gn1.shp

Negeer de waarschuwingen. Die hebben te maken met formaten en tekens die niet ondersteund worden.

Hier kan je alle opties van het ogr2ogr commando vinden. Merk op dat je dit commando kan gebruiken als je tegelijkertijd herprojectie en conversie wilt doen. In dat geval gebruik je -t_srs. Je kan nog veel meer met dit commando, maar dat kunnen we helaas niet behandelen in deze beginnerscursus.

2. Open het bestand in Google Earth om het resultaat te bekijken.

4.3. Ruimtelijke queries met vectordata

Voor onze kaart van Delft willen we de volgende GIS analyses doen:

Selecteer de gemeente Delft van de gemeentenkaart en sla het op in een nieuwe shapefile;
Gebruik de gemeentegrens van Delft om de wegen uit te knippen, zodat we alleen de wegen van Delft hebben.

Met een ruimtelijke querie kunnen we een object uit een vectorlaag selecteren.

1. Welk attribuut wordt in de gemeentenkaart gebruikt voor de gemeentenaam? Gebruik ogrinfo om het antwoord te vinden.

2. Voer het volgende commando uit:

ogr2ogr
-f "ESRI Shapefile" -where GM_NAAM='Delft' 
-a_srs EPSG:28992 delft.shp gem_2011_gn1.shp

Dit commando slaat het object met GM_NAAM Delft op in een bestand met de naam delft.shp. Het argument -a_srs EPSG:28992 gebruiken we om

de Amersfoort/RD Nieuw projectie toe te wijzen aan het uitvoerbestand. Met het argument -f definiëren we het uitvoerformaat.

3. We kunnen de roadsreprojected.shp vector nu knippen, zodat we alleen de wegen in de gemeente Delft over houden. We gebruiken hiervoor het -clipsrc argument:

ogr2ogr -clipsrc delft.shp roadsdelft.shp

roadsreprojected.shp

4. Open nu in QGIS het geherprojecteerde hoogtemodel (dem_rd.tif), de geknipte wegenkaart (roadsdelft.shp) en de gemeentegrens van Delft ( delft.shp).

4.4. Converteer Comma Separated Values (CSV) bestanden

Soms wil je de coordinaten in een tekstbestand converteren. Bijvoordbeeld een ASCII bestand uit een spreadsheet programma. Hier gaan we de coordinaten in het CSV bestand locations.csv herprojecteren en opslaan in een nieuw ASCII bestand met de naam locations_reprojected.csv.

Het is goed gebruik om eerst de inhoud van een CSV bestand in een tekst editor (zoals Kladblok) te bekijken en uit te vinden waar de coordinaten staan en welk teken gebruikt wordt om de kolommen te scheiden. Dat is namelijk niet altijd een komma en hangt soms af van de taalinstellingen die het spreadsheetprogramma gebruikt tijdens het exporteren naar CSV formaat.

1. Bekijk de inhoud van locations.csv . Je kan het type commando gebruiken, zoals je hebt geleerd in de tutorial over de opdrachtregelinterface.

We zien dat de coordinaten in de lat/lon kolommen staan. Dat betekent dat we EPSG:4326 moeten gebruiken. De kolommen worden gescheiden met komma's. In de laatste kolom staat tekst, herkenbaar aan de dubbele aanhalingstekens.

Om de projectie van het CSV bestand te veranderen moeten we eerst een virtuele databron maken d.m.v. een XML bestand.

2. Open Kladblok and typ/kopieer de XML code die hieronder staat. Gebruik drie spaties voor het inspringen.

<OGRVRTDataSource>
   <OGRVRTLayer name="locations">
      <SrcDataSource>locations.csv</SrcDataSource>
      <GeometryType>wkbPoint</GeometryType>
      <LayerSRS>EPSG:4326</LayerSRS>
      <GeometryField encoding="PointFromColumns" x="lon" y="lat"/>
   </OGRVRTLayer>
</OGRVRTDataSource>

3. Sla het bestand op als locations.vrt in de gdal_exercises map.

Laten we het XML bestand nader bekijken:

<OGRVRTLayer name="locations"> moet overeenkomen met <SrcDataSource>locations.csv</SrcDataSource>
<LayerSRS>EPSG:4326</LayerSRS> moet overeenkomen met de EPSG code van de coördinaten in het CSV bestand
<GeometryField encoding="PointFromColumns" x="lon" y="lat"/> geeft aan welke kolommen de coördinaten bevatten.

4. Voer het volgende commando uit:

ogr2ogr
-t_srs EPSG:28992 -f "CSV" locations_reprojected.csv locations.vrt -lco
GEOMETRY=AS_XY

In dit voorbeeld converteren we de lat/lon WGS-84 coordinaten van locations.csv naar Amersfoort/RD Nieuw coordinaten en het resultaat wordt opgeslagen in locations_projected.csv.

5. Gebruik Kladblok om het resultaat te bekijken. Wat is er opgeslagen in de kolommen?

Op dezelfde wijze kunnen we het CSV bestand converteren naar een shapefile.

6. Voer het volgende commando uit:

ogr2ogr
-f "ESRI Shapefile" -t_srs
EPSG:28992 locations.shp locations.vrt

7. Bekijk het resultaat met ogrinfo.

8. Visualiseer de shapefile in QGIS door de lokaties over het hoogtemodel, wegenkaart en gemeentegrens van Delft te plotten. Maak er een mooie kaart van.

9. Converteer de lokaties nu naar het Google KML formaat en bekijk in Google Earth wat de lokaties in Delft zijn.

5. Batch conversie

Desktop GISprogramma's zijn nuttig voor GISbewerkingen, maar minder bruikbaar als we dezelfde taken moeten herhalen voor veel data. In dat geval kan je veel tijd winnen door een script te maken.

We gaan hier de resultaten van een landgebruiksmodel van Dublin gebruiken. Deze data is in het IDRISI raster formaat (.rst) , met voor elk jaar tussen 1990 en 2030 een rasterbestand. Onze taak is om alle rasters te converteren naar GeoTIFF (.tif ) formaat.

1. Unzip landuse.zip in een nieuwe directory.

Voor het converteren van elk afzonderlijk rasterbestand zouden we het volgende commando gebruiken:
gdal_translate -of GTiff 01_State19900101.rst 01_State19900101.tif

We gaan nu een batchbestand (zie tutorial Opdrachtregelinterface) maken met een iteratie zodat we alle bestanden in de directory kunnen converteren.

2. Open een teksteditor, bijvoorbeeld Kladblok

3. Typ/kopieer de volgende code:

for %%f in (*.rst) do (


       echo %%~nf
       gdal_translate -of GTiff %%f %%~nf.tif

)

4. Sla het batchbestand op als rst2tif.bat in de map met de landgebruiksrasters (vergeet niet om de bestandsextensie te veranderen als je Kladblok gebruikt, dat is een klassieke fout!).

Probeer de code te begrijpen. Dit is een zogenaamde "for" loop waarmee we itereren over all *.rst files in de directory. %%f is de variabele met de bestandsnaam voor ieder bestand in de directory. Met echo tonen we iets op het scherm. Hier laten we %%~nf zien, wat het deel van de bestandsnaam voor de punt is. Vervolgens gebruiken we het gdal_translate commando met als uitvoerformaat GeoTIFF. Aan het eind van de regel voegen we .tif toe als bestandsextensie.

5. Voer het batchbestand uit. Typ:
rst2tif <ENTER>

6. Bekijk het resultaat.

6. Conclusie

In deze cursus heb je geleerd om:

Informatie op te vragen uit geografische data
Geografische data te herprojecteren
Rastereigenschappen te wijzigen
Rasterformaten te converteren
Vectorformaten te converteren
Ruimtelijke queries toe te passen op vectoren
Bestanden met door komma's gescheiden waarden (CSV bestanden) te converteren
Batchconversie uit te voeren