Hello everyone, I haven't been writing for quite some time because I've been very busy last months but from now on I'll try to write as often as I can.
Today I'd like to talk a little bit more about GDAL library, WMS service and how to extract different type of images from WMS layers.
Before you keep reading please let me tell you that ever since I started learning more about GIS I've always prefered vector data before raster and I've honestly never used GDAL commands before.
However, now I came up with some project idea where I should find a way to handle some orthophoto imagery and do some raster image processing.
That's why I decided to learn more about it and to share with you some things that I've learned while playing with it.
Let's first give short explanation of terms that will be used here.
GDAL (Geospatial Data Abstraction Library) is a library for reading and writing raster geospatial data formats. It's released by the OSGEO (Open Source Geospatial Foundaiton) under X/MIT type of Open Source license. OGR library is part of the GDAL source tree and it provides similair capability for vector data. You can read more about GDAL/OGR on its official site.
If you want to try by yourself examples described in this post then you should first install GDAL. You are free to choose the way to do it but if you are to lazy to figure out how to do it then this should do the job for you:
tar xvfz gdal-1.9.2.tar.gz
sudo make install
WMS (Web Map Service) is OGC standard for serving georeferenced map images over the Internet. There are many WMS Servers that can serve WMS layers but probably the most often used are Geoserver and Mapserver.
WMS requests support this operations:
GetCapabilities - retrieves metadata about the service (e.g. map image format, WMS version, map bounding box, CORS, ..)
GetMap - retrieves map image for defined area and content
GetFeatureInfo (optional) - retrieves geometry and attribute data for specific location on map
GetLegendGraphic (optional) - returns map legend
To help you get better idea of one simple GetMap request we can try it on some real world example. For this purposes I will use WMS layer which contains DOF (Digital Orthopohoto) imagery of my country (Croatia).
If you just copy and paste this url in your browser you should see small image (256x256 pixels) in left upper corner of the screen. It's response for your GetMap request (if you don't see any image then it's probably because that server is not in function anymore, in that case please write a comment below this post so I can change it).
If you take a closer look at the defined url you can see that we set different parameters there. Here is short description for each of them.
service - service name
request - operation name
version - service version
layers - layers to display
srs - Spatial Reference System for map output
format - format for the map output
width - width of map output (in pixels)
height - height of map output (in pixels)
bbox - bounding box for map extent (minx,miny,maxx,maxy)
Those are just some of parameters you can define for WMS GetMap request but for purpose of this example there is no need for additional parameters.
So, as you can see, we defined set of parameters to get jpeg image of specific area within some WMS layer. Although this can often be very useful and interesting sometimes we need better way to access some WMS service and its data.
For example, fif you need to save response images then you probably don't want to do it manually by typing url in your browser every time and then saving it.
As you can assume it can easily be done with Python.
One way for it is to use urllib library and then directly call a WMS and write the response out to a file. However, I'd like to present you another approach this time which provides you much more possibilites.
Accessing web image services using OWSLib
If you haven't found it already, there is great library called OWSLib.
As it says in its description - "OWSLib is a Python package for client programming with Open Geospatial Consortium (OGC) web service (hence OWS) interface standards, and their related content models."
It can be used for some other OGC services too but in this example we will use it just for working with WMS service.
Here is simple code where you can see how to do previous GetMap request using OWSLib library.
from owslib.wms import WebMapService
wms = WebMapService('http://ganimed.geoportal.dgu.hr/cwms', version='1.1.1')
img = wms.getmap(layers=['DOF'],
srs = 'EPSG:3765',
bbox = (399609, 4887306, 400326, 4888023),
size = (256,256),
format = 'image/jpeg',
transparent = True
out = open('dof_img.jpg', 'wb')
This can be very useful when you develop some web applications and need to make different WMS requests based on your user inputs (e.g. bounding box) and similair.
However, please visit OWSLib official documentation to get better idea of what you can all do with this library.
Accessing web image services using WMS in GDAL
Beside previously described approach, accessing several different types of web image services is also possible using the WMS format in GDAL.
To do that, you first need to create local service description XML file.
For example, if you want to get same image as before (when we used Python) then your XML file should look somehow like this:
Now when you have XML file created you can run GDAL commands on it.
For the beginning let's try to get information about our file:
You should see output with set of different informations (e.g. driver, size, coordinate system, ..)
However, we are probably more interesting in doing something more specific with our file.
Here you can find all utility programs that are distributed with GDAL (http://www.gdal.org/gdal_utilities.html)
Probably the most often used command is gdal_translate. We will use it here to get image output (.jpeg file) of our defined area based on previously created XML file.
gdal_translate -of JPEG -outsize 256 256 gdal_wms_dgu.xml dof_img2.jpg
As you can probably guess this command gives us JPEG image output (256x256 pixels) based on our XML file.
If you compare it to our dgu_dof.jpg image that we got by running our python script then you should see that they are the same.
Please note that you need to explicity define image output because by default it is set to GeoTIFF (GTiff) so instead of .jpg file you would get .tif file.
Actually, let's try to do that because if you are working with map images and GIS applications then you would probably need GeoTIFF images before jpeg.
gdal_translate -of GTiff -outsize 256 256 gdal_wms_dgu.xml dof_geotiff.tif
Now if you want you can get more info about our output image by running:
To find out more possible operations on raster data (like resampling and rescaling pixels, setting ground control points on output images,...) please visit http://www.gdal.org/gdal_translate.html
At the end, I'd just like to mention one more thing about running GDAL commands within Python.
First, it's important to mention that you can use GDAL within Python if you install it as Python Package https://pypi.python.org/pypi/GDAL/.
Honestly as I've just started exploring this raster image processing with GDAL I didn't have time to test how to use the Python GDAL/OGR API but I definitely willl and maybe I'll write new post then.
However, I'd like to remember you that Python has subprocess module which can be used for accessing system commands.
So if you want to run previous GDAL command within python script you can always do it like this:
subprocess.call(["gdal_translate", "-of", "GTiff", "-outsize", "250", "250", "gdal_wms_dgu.xml", "dof_img2.tif"])
That would be all for this post, I hope some of you will find it useful or at least it will help you get an idea for some of you future tasks and projects. I will continue studying this field and I hope I will write more about it.
Feel free to comment and ask anything you want below.