Page 1 of 1 (19 posts)

  • talks about »
  • gdal

Tags

Last update:
Wed Aug 20 07:40:09 2014

A Django site.

QGIS Planet

GDAL/OGR 1.11.0 released

The new version 1.11.0  of GDAL/OGR (http://www.gdal.org/) which offers major new features has been released. GDAL/OGR is a C++ geospatial data access library for raster and vector file formats, databases and web services.  It includes bindings for several languages, and a variety of command line tools.

Highlights:

More complete information on the new features and fixes in the 1.11.0 release can be found at http://trac.osgeo.org/gdal/wiki/Release/1.11.0-News

The new release can be downloaded from:

The post GDAL/OGR 1.11.0 released appeared first on GFOSS Blog | GRASS GIS Courses.

OSGeo Code Sprint, Vienna

This is how OSGeo happens.  These are the folk who bring us a lot of that open-source geo-spatial goodness. You can follow the code sprint on Twitter using the hashtags #csprint and #viennacodesprint14

 


50th ICA-OSGeo Lab established at Fondazione Edmund Mach (FEM)

We are pleased to announce that the 50th ICA-OSGeo Lab has been established at the GIS and Remote Sensing Unit (Piattaforma GIS & Remote Sensing, PGIS), Research and Innovation Centre (CRI), Fondazione Edmund Mach (FEM), Italy. CRI is a multifaceted research organization established in 2008 under the umbrella of FEM, a private research foundation funded by the government of Autonomous Province of Trento. CRI focuses on studies and innovations in the fields of agriculture, nutrition, and environment, with the aim to generate new sharing knowledge and to contribute to economic growth, social development and the overall improvement of quality of life.

The mission of the PGIS unit is to develop and provide multi-scale approaches for the description of 2-, 3- and 4-dimensional biological systems and processes. Core activities of the unit include acquisition, processing and validation of geo-physical, ecological and spatial datasets collected within various research projects and monitoring activities, along with advanced scientific analysis and data management. These studies involve multi-decadal change analysis of various ecological and physical parameters from continental to landscape level using satellite imagery and other climatic layers. The lab focuses on the geostatistical analysis of such information layers, the creation and processing of indicators, and the production of ecological, landscape genetics, eco-epidemiological and physiological models. The team pursues actively the development of innovative methods and their implementation in a GIS framework including the time series analysis of proximal and remote sensing data.

The GIS and Remote Sensing Unit (PGIS) members strongly support the peer reviewed approach of Free and Open Source software development which is perfectly in line with academic research. PGIS contributes extensively to the open source software development in geospatial (main contributors to GRASS GIS), often collaborating with various other developers and researchers around the globe. In the new ICA-OSGeo lab at FEM international PhD students, university students and trainees are present.

PGIS is focused on knowledge dissemination of open source tools through a series of courses designed for specific user requirement (schools, universities, research institutes), blogs, workshops and conferences. Their recent publication in Trends in Ecology and Evolution underlines the need on using Free and Open Source Software (FOSS) for completely open science. Dr. Markus Neteler, who is leading the group since its formation, has two decades of experience in developing and promoting open source GIS software. Being founding member of the Open Source Geospatial Foundation (OSGeo.org, USA), he served on its board of directors from 2006-2011. Luca Delucchi, focal point and responsible person for the new ICA-OSGeo Lab is member of the board of directors of the Associazione Italiana per l’Informazione Geografica Libera (GFOSS.it, the Italian Local Chapter of OSGeo). He contributes to several Free and Open Source software and open data projects as developer and trainer.

Details about the GIS and Remote Sensing Unit at http://gis.cri.fmach.it/

Open Source Geospatial Foundation (OSGeo) is a not-for-profit organisation founded in 2006 whose mission is to support and promote the collaborative development of open source geospatial technologies and data.

International Cartographic Association (ICA) is the world authoritative body for cartography and GIScience. See also the new ICA-OSGeo Labs website.

Processing Landsat 8 data in GRASS GIS 7: Import and visualization

banner_landsat_rgb

For our analysis example, we’ll obtain (freely – thanks, NASA and USGS!) a Landsat 8 scene from http://earthexplorer.usgs.gov/

First of all, you should register.

1. Landsat 8 download procedure

1. Enter Search Criteria:

  • path/row tab, enter Type WRS2: Path: 16, Row: 35
  • Date range: 01/01/2013 – today
  • Click on the “Data sets >>” button

2. Select Your Data Set(s):

  • Expand the entry + Landsat Archive
        [x] L8 OLI/TIRS
  • Click on the “Results >>” button

(We jump over the additional criteria)

4. Search Results

From the resulting list, we pick the data set:

earthexplorer_selection_lsat8Entity ID: LC80160352013134LGN03
Coordinates: 36.04321,-79.28696
Acquisition Date: 14-MAY-13

Using the “Download options”, you can download the data set (requires login). Select the choice:
[x] Level 1 GeoTIFF Data Product (842.4 MB)

You will receive the file “LC80160352013134LGN03.tar.gz”.

2. Unpacking the downloaded Landsat 8 dataset

To unpack the data, run (or use a graphical tool at your choice):

tar xvfz LC80160352013134LGN03.tar.gz

A series of GeoTIFF files will be extracted: LC80160352013134LGN03_B1.TIF, LC80160352013134LGN03_B2.TIF, LC80160352013134LGN03_B3.TIF, LC80160352013134LGN03_B4.TIF, LC80160352013134LGN03_B5.TIF, LC80160352013134LGN03_B6.TIF, LC80160352013134LGN03_B7.TIF, LC80160352013134LGN03_B8.TIF, LC80160352013134LGN03_B9.TIF, LC80160352013134LGN03_B10.TIF, LC80160352013134LGN03_B11.TIF, LC80160352013134LGN03_BQA.TIF

We may check the metadata with “gdalinfo“:

gdalinfo LC80160352013134LGN03_B1.TIF
Driver: GTiff/GeoTIFF
Files: LC80160352013134LGN03_B1.TIF
Size is 7531, 7331
Coordinate System is:
PROJCS["WGS 84 / UTM zone 17N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
...
Pixel Size = (30.000000000000000,-30.000000000000000)
...

3. Want to spatially subset the Landsat scene first?

If you prefer to cut out a smaller area (subregion), check here for gdal_translate usage examples.

4. Import into GRASS GIS 7

Note: While this Landsat 8 scene covers the area of the North Carolina (NC) sample dataset, it is delivered in UTM rather than the NC's state plane metric projection. Hence we preprocess the data first in its original UTM projection prior to the reprojection to NC SPM.

Using the Location Wizard, we can import the dataset easily into a new location (in case you don't have UTM17N not already created earlier):

grass70 -gui

grass7_loc_wizard1
grass7_loc_wizard2
grass7_loc_wizard3
grass7_loc_wizard4
grass7_loc_wizard5
grass7_loc_wizard6
grass7_loc_wizard7
grass7_loc_wizard8
grass7_loc_wizard9

 

 

 

Now start GRASS GIS 7 and you will find the first band already imported (the others will follow shortly!).

For the lazy folks among us, we can also create a new GRASS GIS Location right away from the dataset on command line:

grass70 -c LC80160352013134LGN03_B10.TIF ~/grassdata/utm17n

5. Importing the remaining Landsat 8 bands

The remaining bands can be easily imported with the raster import tool:

grass7_import1

The bands can now be selected easily for import:

grass7_import2

  • Select "Directory" and navigate to the right one
  • The available GeoTIFF files will be shown automatically
  • Select those you want to import
  • You may rename (double-click) the target name for each band
  • Extend the computation region accordingly automatically

Click on "Import" to get the data into the GRASS GIS location. This takes a few minutes. Close the dialog window then.

In the "Map layers" tab you can select the bands to be shown:

grass7_visualize1

6. The bands of Landsat 8

(cited from USGS)

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) images consist of nine spectral bands with a spatial resolution of 30 meters for Bands 1 to 7 and 9. New band 1 (ultra-blue) is useful for coastal and aerosol studies. New band 9 is useful for cirrus cloud detection. The resolution for Band 8 (panchromatic) is 15 meters. Thermal bands 10 and 11 are useful in providing more accurate surface temperatures and are collected at 100 meters. Approximate scene size is 170 km north-south by 183 km east-west (106 mi by 114 mi).

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS)Launched February 11, 2013 Bands Wavelength (micrometers) Resolution (meters)
Band 1 - Coastal aerosol 0.43 - 0.45 30
Band 2 - Blue 0.45 - 0.51 30
Band 3 - Green 0.53 - 0.59 30
Band 4 - Red 0.64 - 0.67 30
Band 5 - Near Infrared (NIR) 0.85 - 0.88 30
Band 6 - SWIR 1 1.57 - 1.65 30
Band 7 - SWIR 2 2.11 - 2.29 30
Band 8 - Panchromatic 0.50 - 0.68 15
Band 9 - Cirrus 1.36 - 1.38 30
Band 10 - Thermal Infrared (TIRS) 1 10.60 - 11.19 100
Band 11 - Thermal Infrared (TIRS) 2 11.50 - 12.51 100

* TIRS bands are acquired at 100 meter resolution, but are resampled to 30 meter in delivered data product.

7. Natural color view (RGB composite)

Due to the introduction of a new "Cirrus" band (#1), the BGR bands are now 2, 3, and 4, respectively. See also "Common band combinations in RGB" for Landsat 7 or Landsat 5, and Landsat 8.

From Digital Numer (DN) to reflectance:
Before creating an RGB composite, it is important to convert the digital number data (DN) to reflectance (or optionally radiance). Otherwise the colors of a "natural" RGB composite do not look convincing but rather hazy (see background in the next screenshot). This conversion is easily done using the metadata file which is included in the data set with i.landsat.toar:

grass7_landsat_toar0
grass7_landsat_toar1
grass7_landsat_toar2
grass7_landsat_toar3

Now we are ready to create a nice RGB composite:

grass7_landsat_rgb0

grass7_landsat_rgb1

Select the bands to be visually combined:

grass7_visualize2

... and voilà!

grass7_landsat_rgb2

8. Applying the Landsat 8 Quality Assessment (QA) Band

One of the bands of a Landsat 8 scene is named "BQA" which contains for each pixel a decimal value representing a bit-packed combination of surface, atmosphere, and sensor conditions found during the overpass.  It can be used to judge the overall usefulness of a given pixel.

We can use this information to easily eliminate e.g. cloud contaminated pixels. In short, the QA concept is (cited here from the USGS page):

Cited from http://landsat.usgs.gov/L8QualityAssessmentBand.php‎

For the single bits (0, 1, 2, and 3):
0 = No, this condition does not exist
1 = Yes, this condition exists.

The double bits (4-5, 6-7, 8-9, 10-11, 12-13, and 14-15) represent levels of confidence that a condition exists:
00 = Algorithm did not determine the status of this condition
01 = Algorithm has low confidence that this condition exists (0-33 percent confidence)
10 = Algorithm has medium confidence that this condition exists (34-66 percent confidence)
11 = Algorithm has high confidence that this condition exists (67-100 percent confidence).

Detailed bit patterns (d: double bits; s: single bits):
d  Bit 15 = 0 = cloudy
d  Bit 14 = 0 = cloudy
d  Bit 13 = 0 = not a cirrus cloud
d  Bit 12 = 0 = not a cirrus cloud
d  Bit 11 = 0 = not snow/ice
d  Bit 10 = 0 = not snow/ice
d  Bit 9 = 0 = not populated
d  Bit 8 = 0 = not populated
d  Bit 7 = 0 = not populated
d  Bit 6 = 0 = not populated
d  Bit 5 = 0 = not water
d  Bit 4 = 0 = not water
s  Bit 3 = 0 = not populated
s  Bit 2 = 0 = not terrain occluded
s  Bit 1 = 0 = not a dropped frame
s  Bit 0 = 0 = not fill

Usage example 1:  Creating a mask from a bitpattern

We can create a cloud mask (bit 15+14 are set) from this pattern:
cloud:   1100000000000000

Using the Python shell tab, we can easily convert this into the corresponding decimal number for r.mapcalc:

Cited from http://landsat.usgs.gov/L8QualityAssessmentBand.php‎

Welcome to wxGUI Interactive Python Shell 0.9.8

Type "help(grass)" for more GRASS scripting related information.
Type "AddLayer()" to add raster or vector to the layer tree.

Python 2.7.5 (default, Aug 22 2013, 09:31:58) 
[GCC 4.8.1 20130603 (Red Hat 4.8.1-1)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> int(0b1100000000000000)
49152

Using this decimal value of 49152, we can create a cloud mask:

# set NULL for cloudy pixels, 1 elsewhere:
r.mapcalc "cloudmask = if(LC80160352013134LGN03_BQA == 49152, null(), 1 )"

# apply this mask
r.mask cloudmask

In our sample scene, there are only tiny clouds in the north-east, so no much to be seen. Some spurious cloud pixels are scattered over the scene, too, which could be eliminated (in case of false positives) or kept.

Usage example 2:  Querying the Landsat 8 BQA map and retrieve the bitpattern

Perhaps you prefer to query the BQA map itself (overlay the previously created RGB composite and query the BSA map by selecting it in the Layer Manager). In our example, we query the BQA value of the cloud:

Using again the Python shell tab, we can easily convert the decimal number (used for r.mapcalc) into the corresponding binary representation to verify with the  table values above.

>>> x=61440
>>> print(bin(x & 0xffffffff))
0b1111000000000000

Hence, bits 15,14,13, and 12 are set: cloudy and not a cirrus cloud. Looking at the RGB composite we tend to agree :-) Time to mask out the cloud!

wxGUI menu >> Raster >> Mask [r.mask]

Or use the command line, as shown already above:

# remove existing mask (if active)
r.mask -r

# set NULL for cloudy pixels, 1 elsewhere:
r.mapcalc "cloudmask = if(LC80160352013134LGN03_BQA == 61440, null(), 1 )" --o

# apply the new mask
r.mask cloudmask

The visual effect in the RGB composite is minimal since the cloud is white anyway (as NULL cells, too). However, it is relevant for real calculations such as NDVI (vegetation index) or thermal maps.

We observe dark pixels around the cloud originating from thin clouds. In a subsequent identification/mask step we may eliminate also those pixels with a subsequent filter.

Installing python GDAL into a python virtualenv in OSX

Yes I know its a bit wierd to see the word ‘OSX’ in an article by me…I recently bought a Mac for testing InaSAFE and fooling around with FOSSGIS software on OSX. Besides the 2 or three days it took to get it set up as a development machine for FOSSGIS stuff (as compared to... Read more »

Installing python GDAL into a python virtualenv in OSX

Yes I know its a bit wierd to see the word 'OSX' in an article by me...I recently bought a Mac for testing InaSAFE and fooling around with FOSSGIS software on OSX. Besides the 2 or three days it took to get it set up as a development machine for FOSSGIS stuff (as compared to a few hours - limited mainly by internet speed - on Linux) it hasnt been too bad an experience.

I used William Kyngesburye's excellent packages for installing all the FOSSGIS dependencies (barring QGIS which of course I wanted to build myself). I've been trying to see if I can use OSX as a usable platform for doing GeoDjango work and ran into a few problems:

  • GDAL python package needs compilation when used in a python virtualenv (which of course I use)
  • William's packages provide PostGIS 2.0 whereas GeoDjango really only works out of the box with 1.5 still

In this post I am outlining how I got GDAL to build in my venv. Note that I assume you have Williams frameworks and xcode / compilation toolchain in place.

source venv/bin/activate
pip install --no-install GDAL
cd venv/build/GDAL
python setup.py build_ext \
  --gdal-config=/usr/local/gdal/1.9/bin/gdal-config \
  --library-dirs=/Library/Frameworks/GDAL.framework/Versions/1.9/unix/lib/ \
  --include-dirs=/Library/Frameworks/GDAL.framework/Versions/1.9/Headers/
pip install --no-download GDAL

Now verify it is installed:

pip freeze

You should see GDAL listed. Another check:

>>> from osgeo import gdal
>>> gdal.__version__
'1.9.2'

Confirm it is really using your module in your venv:

>>> gdal.__file__
'/Users/timlinux/dev/python/catalogue/venv/lib/python2.7/site-packages/osgeo/gdal.pyc'

WFS to PostGIS in 1 Step

This is an update to my previous post “WFS to PostGIS in 3 Steps”. Thanks to Even Rouault’s comments and improvements to GDAL, it is now possible to load Latin1-encoded WFS (like the one by data.wien.gv.at) into PostGIS in just one simple step.

To use the following instructions, you’ll have to get the latest GDAL (release-1600-gdal-mapserver.zip)

You only need to run SDKShell.bat to set up the environment and ogr2ogr is ready for action:

C:\Users\Anita>cd C:\release-1600-gdal-mapserver
C:\release-1600-gdal-mapserver>SDKShell.bat
C:\release-1600-gdal-mapserver>ogr2ogr -overwrite -f PostgreSQL PG:"user=myuser password=mypassword dbname=wien_ogd" "WFS:http://data.wien.gv.at/daten/geoserver/ows?service=WFS&request=GetFeature&version=1.1.0&typeName=ogdwien:BEZIRKSGRENZEOGD&srsName=EPSG:4326"

Thanks everyone for your comments and help!


Creating coloured rasters with GDAL

One of the most popular posts on my blog ever is my article about creating GDAL hillshades: ‘A workflow for creating beautiful relief shaded DEMS using GDAL‘. In the techniques outlined in the aforementioned article, colours are assigned across a continuum to and superimposed on top of a hillshaded model. The GDAL tools in QGIS [...]

Generating contours using GDAL ( via shell or QGIS)

I tell you. It always amazes me how much cool stuff you can do with great open source GIS software these days. GDAL is one of those great open source projects that I have just found a great use for (apart from just opening every raster type under the sun in QGIS).

GDAL has the ability to generate contours from a DEM, something that I have always wanted to try for my town but have never been able due to lack of a good DEMs.

Recently we purchased a set of DEMs that cover a large area as part of a study. Each DEM uses a grid size of 1mx1m. Prefect for generating contours.

Using the GdalTools QGIS plugin.

First make sure that you have the latest version of the GdalTools plugin installed (GdalTools should be installed by default with QGIS. If it’s not, search “Gdal” in the plugin installer). Enable the plugin once it’s installed.

Load the DEM into QGIS using the Load Raster icon.

DEM loaded in QGIS

Now head up to the menu Raster->Extraction->Contour

Raster menu in QGIS

Select the settings that you need. For this DEM I am going to generate 250mm contours.

Contour dialog.

Take note of the text area at the bottom of the dialog as that is the exact command sent to GDAL in order to generate the contours. If you take a copy of that you can run it on the command line for batch processing later.

Hit ok.

250mm contours from the DEM

BAM! :)

Using the command line/shell.

Using QGIS for a one off DEM is fine and dandy but what if you have 3000 DEMs that you need to process. To hell with doing that by hand!

Remember the contour tool in QGIS told us the exact command line args to use, so creating a shell script for automating the process is pretty easy.

for f in *.asc
do
  echo "Processing $f"
 gdal_contour -a ELEV -i 0.25 $f $f-250mm.shp
done

The above code will loop though the current directory and process all the DEMs generating 250mm contours for each one. It saves each new contour file as {filename}-250mm.shp. You will need to change *.asc to whatever format your DEM is in.

Copy the above code into a file somewhere and call it generate_contours.sh. This can then be called from the command line using

sh generate_contours.sh

Running sh on Windows

If you’re a windows user you will need to run OsGeo4W Shell in order to use sh.

Loading OsGeo4w shell.

Once the shell is loaded you can just call:

sh generate_contours.sh

Output from running generate_contours.sh

Final remarks

I ran the above sh script on a folder with about 2500 DEMs and it took around 4 hours to complete the whole folder. Of course performance will vary but it seems pretty quick considering.

Once again the possibly to use open source GIS tools to get my work done is amazing.  No expensive software here.

So far I’m not aware of any line smoothing/generalizing abilities using GDAL/OGR.  Although you can import the contours into GRASS and use that: http://grass.osgeo.org/wiki/V.generalize_tutorial

You can also generate the contours using GDAL via Python but that is a topic for another day.


Filed under: Open Source, qgis Tagged: contours, DEM, FOSSGIS, gdal, gis, grid, ogr, Open Source, osgeo, qgis, Quantum GIS, raster

Selecting GDAL Drivers on the fly with QGIS

There is quite an old ticket in the QGIS bug tracker relating to the need to be able to select which GDAL drivers to use. The issues is this: GDAL in some cases provides multiple drivers for the same image type. For example, JPEG2000 datasets can be opened using ECW's proprietary driver, KAKADU (gotta love that name!), OpenJPEG etc. The problem is that there isn't actually any way to choose which driver GDAL should use for a given image format - it works on a first come, first served bases. This is a bit of a problem sometimes when you actually want to be using a different driver. In some work I did recently for SANSA (the South African Space Agency), I built a custom version of QGIS for them. They are particularly keen on the JPEG2000 format since it is an open standard and provides good compression. However the custom build of QGIS I made for them includes the ECW JP2000 driver which unfortunately chokes and dies (no fault of GDAL, it happens down in the ECW code itself) when trying to open their JP2 imagery. Also we would prefer to use a FOSS driver for JP2 in keeping with the general ethos of being a FOSS project. So all this prompted me to implement the needed support to let you disable GDAL drivers at runtime in QGIS which is illustrated a little by the screenshots below.

 

By default GDAL want to use the ECW JPG2000 driver...(click to enlarge)

By default GDAL want to use the ECW JPG2000 driver...(click to enlarge)

Using the options dialog, you can now disable a driver....(click to enlarge)

Using the options dialog, you can now disable a driver....(click to enlarge)

After the driver is disabled (and restarting QGIS), the alternate JP2 driver is available (click to enlarge)

After the driver is disabled (and restarting QGIS), the alternate JP2 driver is available (click to enlarge)

pixelstats trackingpixel

Improvements to raster performance in QGIS master

Today I committed changes to QGIS master that I started working on at the Lisbon hackfest and finalised this week (along with a unit test yay!). The changes implement improvements to the way in which statistics are gathered for rasters. First a little history: the raster stats code in QGIS is some of the earliest work I have done in QGIS. The code was also hugely inefficient - basically it completely scans every pixel in the file twice! It does this because it extracts minima, maxima, stddev etc. from each band in the raster - and the algorithm used to calculate the stddev requires first calculating the sum of means, hence the second pass. Of course we use the underlying GDAL library to traverse the raster, but we completly bypassed the GDAL stats gathering functions. Now from my testing (see below) the difference of gathering stats in QGIS as opposed to directly in GDAL was not that great. However, the real killer is the fact that QGIS had no reasonable implementation to cache these statistics - either within a session or between sessions. This resulted in quite a lot of  complaints about the slow loading of rasters in QGIS and the spontaneous recollection of stats when making various changes to raster symbology etc.

With the release of 1.7, Radim Blazek implemented 'proper' support for raster providers so that they follow a similar model to the vector providers in QGIS. This opens up room for provider specific optitmisations and for implementing new providers in the future e.g. for PostGIS Raster, Terralib etc. At the Lisbon hackfest, I refactored the provider interface so that it provides a default implementation of stats gathering which can be overloaded by specific providers, each then being able to use best practice for that datasource when gathering stats. I then updated the GDAL Provider to use gdal native calls to fetch statistics.

When GDAl gathers stats, it stores them in a file named after the source raster but with a .aux.xml extension. As a bonus, GDAL also stores histogram data in this file. Thanks to some excellent tips from Etienne (Sky?) on the QGIS mailing list my implementation was further improved so that it actually uses these cached stats from GDAL properly. So let's see how performance is improved with these changes. I tested using a ~650mb tif hillshaded DEM from one of my clients.

 

Raster Load Performance in QGIS

 Added bonus: My changes also use a more efficient approach for obtaining min/max for a band (before all stats are gathered) and also histgram calculations are cached now. So after the initial gathering of stats and histograms, all activities when working with GDAL rasters should be near instant now. By the way the 1 in the graph above is rounded up - in my usage it was near instant in most cases.

I need to still evaluate whether it will be good to backport these changes to the 1.7 stable tree, so for now if you want to benefit from these changes, you need to be using a nightly build or build from source.

Every 6 months we hold a hackfest / QGIS meeting and many of our users generously donate towards it. I hope the above article gives those who do support our meetings a sense of the tangible benefits that arise from these meetings - there have been many other improvements coming into master that arose from the last hackfest that we all get to enjoy! If you would like to help towards the costs of the next hackfest, visit http://qgis.org and click the 'Donate!' button on the left hand side. Have fun!

 

pixelstats trackingpixel

Improvements to raster performance in QGIS master

Today I committed changes to QGIS master that I started working on at the Lisbon hackfest and finalised this week (along with a unit test yay!). The changes implement improvements to the way in which statistics are gathered for rasters. First a little history: the raster stats code in QGIS is some of the earliest work I have done in QGIS. The code was also hugely inefficient - basically it completely scans every pixel in the file twice! It does this because it extracts minima, maxima, stddev etc. from each band in the raster - and the algorithm used to calculate the stddev requires first calculating the sum of means, hence the second pass. Of course we use the underlying GDAL library to traverse the raster, but we completly bypassed the GDAL stats gathering functions. Now from my testing (see below) the difference of gathering stats in QGIS as opposed to directly in GDAL was not that great. However, the real killer is the fact that QGIS had no reasonable implementation to cache these statistics - either within a session or between sessions. This resulted in quite a lot of  complaints about the slow loading of rasters in QGIS and the spontaneous recollection of stats when making various changes to raster symbology etc.

With the release of 1.7, Radim Blazek implemented 'proper' support for raster providers so that they follow a similar model to the vector providers in QGIS. This opens up room for provider specific optitmisations and for implementing new providers in the future e.g. for PostGIS Raster, Terralib etc. At the Lisbon hackfest, I refactored the provider interface so that it provides a default implementation of stats gathering which can be overloaded by specific providers, each then being able to use best practice for that datasource when gathering stats. I then updated the GDAL Provider to use gdal native calls to fetch statistics.

When GDAl gathers stats, it stores them in a file named after the source raster but with a .aux.xml extension. As a bonus, GDAL also stores histogram data in this file. Thanks to some excellent tips from Etienne (Sky?) on the QGIS mailing list my implementation was further improved so that it actually uses these cached stats from GDAL properly. So let's see how performance is improved with these changes. I tested using a ~650mb tif hillshaded DEM from one of my clients.

Raster Load Performance in QGIS

 Added bonus: My changes also use a more efficient approach for obtaining min/max for a band (before all stats are gathered) and also histgram calculations are cached now. So after the initial gathering of stats and histograms, all activities when working with GDAL rasters should be near instant now. By the way the 1 in the graph above is rounded up - in my usage it was near instant in most cases.

I need to still evaluate whether it will be good to backport these changes to the 1.7 stable tree, so for now if you want to benefit from these changes, you need to be using a nightly build or build from source.

Every 6 months we hold a hackfest / QGIS meeting and many of our users generously donate towards it. I hope the above article gives those who do support our meetings a sense of the tangible benefits that arise from these meetings - there have been many other improvements coming into master that arose from the last hackfest that we all get to enjoy! If you would like to help towards the costs of the next hackfest, visit http://qgis.org and click the 'Donate!' button on the left hand side. Have fun!

Another cool open source project – OSGeo-Live

Another cool open source project that I have become a part of (as a QGIS packager and tester) is the OSGeo-Live project.  The OSGeo-Live project is a live DVD/USB/Virtual Machine built on xUbuntu(striped down Ubuntu linux) that has a lot of cool open source geo spatial programs all set up and ready to use.

The OSGeo-Live project contains:

  • Browser clients
  • A small sample of crisis management software
  • All the popular database engines (PostGIS, SpaitalLite etc)
  • Pretty much all the open source desktop GIS apps (QGIS, uDig etc)
  • Open Source GPS navigation apps and globes.
  • A collection of handy spatial tools
  • A ready to go web services ready to try in your browser or desktop GIS.
  • Some sample data to get started with for each project
  • And quick starts for each program.

The full list of software contained on the OSGeo-Live project can be found at http://live.osgeo.org/en/overview/overview.html

This is a good project if you want to get into the OSGeo tools and experiment but don’t want to install them on main machine until you know what you need.

As it is a live DVD/USB/Virtual Machine some apps will run slower than what they do on a native install but overall the speed is usable and good enough for testing.

Even better is that it was born in Australia :)

The project is also commercial supported by a Australian company http://lisasoft.com

So give it a try if you are interested in the OSGeo movement, which you should be if you are reading my blog :)


Filed under: Open Source, qgis Tagged: gdal, gis, mapping, Open Source, osgeo, osgeo-live, OSS, Quantum GIS

Frank Warmerdam to become a Google guy

If you follow the osgeo planet blog aggregator, you may have noticed Frank Warmerdam's recent blog post mentioning that he is off to work for Google. If I think back to the blog posts I have written, GDAL features in a great many of them - its a mainstay tool for me and barely a day goes buy that I don't gdal_translate some file or what-have-you. Hell, I just finished running a batch job that ran for about 3 weeks comprised of nothing more than bash script and GDAL commands.

In QGIS too, GDAL is one of the key underpinnings of our project - it is only now 8 years or so into the project that we are starting to abstract the use of GDAL for reading rasters away a little to make space of other raster provider types. And of course OGR is one of the main ways that our users interact with vector data from within QGIS.

Personally, Frank has also been extremely helpful whenever I have wondered into the #gdal IRC channel looking for help. Reading his blog post made me think that Google doesn't quite realise yet just how lucky they are to get him. If Google has any heart, they will make Frank a 'minister without portfolio' and just tell him to keep working on GDAL and enjoy the free food in their cafetarias. Here's hoping at least that he still gets to spend a little time on the software that so many of use love using and depend on.

pixelstats trackingpixel

3D DXF from MapInfo Tab or ESRI Shape (or anything) using OGR

Note:The following post requires gdal 1.8

A lot of times we need to send/use 3D dxf files, we used to use FME however FME is not free or cheap. So I went looking for a free solution.

If you have gdal 1.8 it’s just one simple command line run using, what is becoming my favorite GIS tool, ogr2ogr.

All you have to do is run:

ogr2ogr -f "DXF" {outFile} {inFile} -zfield {ColumnWithZValue}

So in my case I ran:

ogr2ogr -f "DXF" C:\Temp\contourswarwick.dxf C:\Temp\Contours.TAB -zfield Height

I haven’t fully tested it but {outfile} can be any file ogr supports.

and the output:

Top view of contours

Top view of contour layer

Contours side view

The side view of the above image.

In the words of, the not so great, Charlie Sheen. WINNING!

Happy mapping :D


Filed under: MapInfo, Open Source Tagged: 3D, ESRI, gdal, gis, mapinfo, mapping, ogr, Open Source, OSS

Importing a DBF containing X-Y Values into QGIS

Suppose you have a DBF (.dbf) file containing X and Y values that you want to import and save as a spatial layer.

QGIS doesn’t support direct import of a DBF file as a map layer, however, we can use some command line magic to convert it to a CSV file and then use the Delimited Text plugin to get the job done.

Your DBF file should have an id for each record and fields containing X and Y values. If it has additional fields that should be OK as well.

First convert the DBF to a comma delimited file using ogr2ogr:
ogr2ogr -f CSV my_csv my_data.dbf

If you don’t have ogr2ogr see http://www.gdal.org/index.html.

This will create my_data.csv in the my_csv directory. You are now ready to bring it into QGIS.

Here are the steps to import the CSV:

  1. Start QGIS
  2. If not already enabled, use the plugin manager to enable the Delimited Text plugin
  3. Click on the Delimited Text icon in the Plugin toolbar or choose it from the Plugins menu
  4. Browse to the location of your CSV file
  5. Enter a name for the layer
  6. Under Selected delimiters, check Comma
  7. If your X and Y fields aren’t automatically determined, set them using the drop-down boxes
  8. The sample text should show how the file is being parsed—if it looks right click OK, otherwise adjust the settings
  9. The layer is added to QGIS

Delmited Text plugin ready to import CSV

At this point the layer behaves pretty much like any other QGIS layer. To save it as a shapefile, right click on its name in the legend and choose Save as…

GDAL/OGR 1.8.0 Released

The 1.8.0 release is a major new feature release with the following highlights:

  • New GDAL drivers: GTX, HF2, JPEGLS, JP2OpenJPEG, JPIPKAK, KMLSUPEROVERLAY, LOS/LAS, MG4Lidar, NTv2, OZI, PDF, RASDAMAN, XYZ
  • New OGR drivers: AeronavFAA, ArcObjects , GPSBabel, HTF, LIBKML, MSSQLSpatial, NAS, OpenAir, PDS, PGDump, SOSI, SUA, WFS
  • Significantly improved OGR drivers: DXF, GML
  • New implemented RFCs
  • New utility: gdallocationinfo

More complete information on the new features and fixes can be found in the GDAL/OGR 1.8.0 Release News.


Batch convert a directory of tiffs to ecw

Today I wanted to batch convert a directory of .tiff images to .ecw (MrSid wavelet compressed). Our server has 8 cores so it would be nice to use them all right?

Here is the quick & dirty way I do this kind of job in parallel.

#!/bin/bash
mkdir ecw
for FILE in *.tif
do
  BASENAME=$(basename $FILE .tif)
  OUTFILE=ecw/${BASENAME}.ecw
  echo "Processing: ${BASENAME}.tif"
  if [ -f $OUTFILE ] #skip if exists
  then
    echo "Skipping: $OUTFILE"
  else
    /usr/local/bin/gdal_translate -of ECW -co LARGE_OK=YES $FILE $OUTFILE
  fi
done

The script is extremely simple and is set up so that you can run it multiple times without problems because if looks to see if the output file already exists before trying to write it. If it does exist, it skips straight on to the next image.

To run 8 parallel processes I simply do this at the command prompt (I did mine in a screen session):

./toecw &
./toecw &
./toecw &
./toecw &
./toecw &
./toecw &
./toecw &
./toecw &

Afterwards you can fire up top and watch 'em go!

top - 18:21:04 up  6:41,  4 users,  load average: 10.29, 9.83, 6.69
Tasks: 216 total,   1 running, 215 sleeping,   0 stopped,   0 zombie
Cpu0  : 56.5%us, 22.5%sy,  0.0%ni, 15.7%id,  4.9%wa,  0.0%hi,  0.3%si,  0.0%st
Cpu1  : 53.3%us, 31.6%sy,  0.0%ni,  8.9%id,  6.2%wa,  0.0%hi,  0.0%si,  0.0%st
Cpu2  : 50.7%us, 37.5%sy,  0.0%ni,  4.9%id,  6.6%wa,  0.0%hi,  0.3%si,  0.0%st
Cpu3  : 46.6%us, 38.4%sy,  0.0%ni,  4.9%id,  9.8%wa,  0.0%hi,  0.3%si,  0.0%st
Cpu4  : 44.0%us, 29.8%sy,  0.0%ni,  8.7%id, 17.2%wa,  0.0%hi,  0.3%si,  0.0%st
Cpu5  : 30.7%us, 57.4%sy,  0.0%ni,  1.7%id,  9.6%wa,  0.0%hi,  0.7%si,  0.0%st
Cpu6  : 58.3%us, 23.8%sy,  0.0%ni,  9.4%id,  8.5%wa,  0.0%hi,  0.0%si,  0.0%st
Cpu7  : 46.1%us, 38.6%sy,  0.0%ni, 10.1%id,  4.6%wa,  0.0%hi,  0.7%si,  0.0%st
Mem:  16227956k total, 16144508k used,    83448k free,  1739140k buffers
Swap: 62492832k total,        0k used, 62492832k free, 13383020k cached

  PID USER      PR  NI  VIRT  RES  SHR S %CPU %MEM    TIME+  COMMAND
12717 timlinux  18  -2  197m  85m 5384 D  104  0.5   0:55.49 gdal_translate
12536 timlinux  18  -2  171m  77m 5384 S  102  0.5   1:08.95 gdal_translate
12705 timlinux  18  -2  195m  65m 5384 D  100  0.4   0:52.58 gdal_translate
12737 timlinux  18  -2  194m  64m 5384 D   97  0.4   0:40.78 gdal_translate
12549 timlinux  18  -2  195m 103m 5384 S   95  0.7   1:12.68 gdal_translate
12751 timlinux  18  -2  165m  66m 5384 S   88  0.4   0:37.46 gdal_translate
12561 timlinux  18  -2  166m  67m 5384 D   69  0.4   1:03.91 gdal_translate
12528 timlinux  18  -2  164m  65m 5384 S   16  0.4   0:18.24 gdal_translate

One thing to note - I ran this with the data sitting on a storage array - if your data all lives on a single drive you may have serious IO issues doing the above....

Extracting bands from hdf4_eos to tif

Here is a quick and dirty bash script I wrote to extract all bands from a modis image and save each to a tiff file. The saved bands are named after the sub layer in the original hdf4 file. The script actually generates a second script called batch.sh so you can peek in there to see the undelying gdal_translate commands that are generated.

One day I need to get around to letting QGIS natively support HDF - maybe a good topic for those attending the OSGEO hackfest in NY?

#!/bin/bash
#Set field separator to linefeeds rather than spaces
echo "#!/bin/bash" > batch.sh
chmod +x batch.sh
IFS=$'\n'
for FILE in *.hdf
do
  echo $FILE
  LAYERS=$(gdalinfo $FILE | grep SUBDATASET | grep NAME  | sed 's/[A-Z0-9\_]*=//g')
  for LAYER in $LAYERS
  do
    LAYER=$(echo $LAYER | sed 's/  //g')
    NEWFILE=$(echo $LAYER | \
    sed 's/  //g' | sed 's/"//g' | sed 's/ //g'| sed 's/:/./g' ).tif
    CMD="gdal_translate -of GTiff '${LAYER}' $NEWFILE"
    echo $CMD >> batch.sh
  done
done

exec ./batch.sh

  • Page 1 of 1 ( 19 posts )
  • gdal

Back to Top

Sponsors