First off, if you haven’t yet, get the basic packages for compiling software, etc. (I enter my username into the sudoers file, similar to Ubuntu, so that sudo works, without being ‘root’ all the time)

~$ sudo yum install gcc gcc-c++ gcc-gfortran flex byacc cmake ccmake \

           ncurses-devel bison wget gsl gsl-devel expat-devel

Next add the available geospatial packages, as well as several python and Qt libraries that we’ll need for GRASS and QGIS packages that we can install without compiling:

~$sudo yum install proj proj-epsg proj-devel geos libgeotiff libgeotiff-devel \
     wxGTK wxGTK-devel swig python-psycopg2 numpy \
     PyQt4 PyQt4-devel sip sip-devel wxPython wxPython-devel

My goal is to have available in GRASS and QGIS both Spatialite based GIS layers and ECW compressed rasters. These libraries are not available on Fedora, so I’ll need to add them by hand and then compile GDAL to include support for these formats. Next I’ll need to work my way up the toolchain, compiling both GRASS and QGIS to use this new GDAL. We begin with spatialite. Get the source:

~$ wget http://www.gaia-gis.it/spatialite-2.3.1/libspatialite-amalgamation-2.4.0.tar.gz
~$ tar xvzf libspatialite-amalgamation-2.4.0.tar.gz
~$ cd libspatialite-amalgamation-2.4.0

At this point we already bump into a problem with the lack of static binaries in Fedora packaging. Spatialite requires the  geos and proj packages. But it looks for *.la files which are missing from all Fedora libraries now. It seems that libtool reads these *.la text files to find what libraries should be linked. In the absence of a suitable *.la file, when we try to compile spatialite, libtool causes a failure. So here’s the workaround I used. I just create three text files with the definitions needed, and then spatialite happily compiles. In /usr/lib you’ll need libgeos_c.la, libgeos.la, and libproj.la. Thanks to Sandro Furieri on the SpatiaLite Users Forum for this tip.  They should look like this:

~$ sudo vi /usr/lib/libgeos_c.la
          dlname='libgeos_c.so'
          library_names='libgeos_c.so'
          inherited_linker_flags=''
          dependency_libs=' -L/usr/lib /usr/lib/libgeos.la'
          weak_library_names=''
          current=7
          age=6
          revision=2
          installed=yes
          shouldnotlink=no
          dlopen=''
          dlpreopen=''
          libdir='/usr/lib'

For each of the additional *.la files change the ‘dlname=’ and ‘library_names=’  entries as required.
Now go ahead and

~$./configure
~$make && sudo make install

.
Next step is the ECW libraries.  This has been faily well documented already. First you’ll have to register on the ERDAS site to download their SDK, then compile and install it. What you will need, after compiling the ECW library, is to add a *.conf file in /etc/ld.so.conf.d so that ldconfig will “find” the new library.

~$sudo vi /etc/ld.so.conf.d/gdal.conf

       /usr/local/lib
~$sudo ldconfig

Now we’re ready for GDAL, etc.  First download the sources:

~$wget http://download.osgeo.org/gdal/gdal-1.7.2.tar.gz
~$wget http://grass.osgeo.org/grass64/source/grass-6.4.0RC6.tar.gz

(or a more recent source from svn)

~$wget http://download.osgeo.org/qgis/src/qgis_1.5.0.tar.gz

Unzip each tarball and, beginning with GDAL we trip over the next change in Fedora. Along with the “no static binaries” policy, the linker ‘ld’ no longer links libraries that are not explicitly listed in the gcc compile statement. Previously, if one binary linked to a second that linked to a third, then there was “implicit” linking that brought in the third binary. This is no longer the case. The error that pops up with gdal concerns the pthread library. GDAL doesn’t use this inherently, but it is used somewhere along the compile chain. Again as a “quick and dirty” workaround I stuck in statements to have GDAL explicitly link to pthread as follows:

~$ cd gdal-1.7.2
~$LIBS="-lpthread" LDFLAGS="-lpthread" ./configure --with-geos=yes --with-sqlite3 --with-ecw=/usr/local --with-python --with-spatialite=/usr/local
~$make
~$sudo make install

With this done, we’re ready for a straightforward compile of GRASS and QGIS:

~$ cd ../grass-6.4.svn_src_snapshot_2010_08_07 

 (or whatever svn snapshot you downloaded…) 

./configure --prefix=/usr/local/grass-6.4 \
     --with-proj-share=/usr/share/proj \
     --with-postgres --enable-largefile  --with-sqlite \
     --with-freetype  --with-freetype-includes=/usr/include/freetype2 \
     --with-python --with-wxwidgets --with-cxx
make && sudo make install
sudo ln -s /usr/local/grass-6.4/bin/grass64 /usr/bin/

You most likely will also want the gdal-grass package:

~$wget http://download.osgeo.org/gdal/gdal-grass-1.4.3.tar.gz
~$/configure --with-grass=/usr/local/grass-6.4/grass-6.4.0svn --with-gdal=/usr/local/bin/gdal-config
~#make && sudo make install

Before moving on to QGIS, we’ll get the R statistics package setup.

~$sudo yum install R-core rpy R-devel xml2 libxml2 libxml2-devel

I’ll now add the spatial packages to R. I like to keep all R packages under the system directory /usr/lib/R. However this directory is read only. Since this is my machine, I allow myself to change permissions, and then I can install additional packages without creating a separate R package directory in my home dir.

~$sudo chgrp -R <my username> /usr/lib/R
~$sudo chmod -R g+w /usr/lib/R

Now open an R session and run:

~$R
> install.packages('sp', dependencies='true')
> install.packages('spgrass6', dependencies='true')
>q()

Next, for the rgdal package. R needs to locate the grass libraries. The R rpm packages which I installed above assume that the grass libs are in /usr/local/lib. But I installed GRASS myself elsewhere, so a quick soft link solves that issue:

~$sudo ln -s /usr/local/grass-6.4/grass-6.4.0svn/lib/libgrass_* /usr/local/lib
~$wget http://cran.r-project.org/src/contrib/rgdal_0.6-28.tar.gz
~$R CMD INSTALL rgdal_0.6-28.tar.gz -l /usr/lib/R/library/

and that’s done.

And finally QGIS:

~$ cd ../qgis-1.5.0
~$mkdir build; cd build
~$cmake -DCMAKE_INSTALL_PREFIX:PATH="/usr/local/qgis-1.5.0"  -DGRASS_PREFIX:PATH="/usr/local/grass-6.4/grass-6.4.0svn" -L ..

(Take note of the ‘..’ at the end of the cmake line). The -L option to cmake outputs a list of the config options. Scan thru it to make sure everything is located on your system. Then:

make && sudo make install
sudo ln -s /usr/local/qgis-1.5.0/bin/qgis /usr/bin

At this point you should be back in Kansas

If you find this learning material worthwhile, or have any comments for improvement, please let me know. I expect that I’ll continue polishing up the material, adding new presentations, and perhaps additional data sources. So check back.

Licensing details:

GIS Training with QGIS by Micha Silver is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

All GPS manufacturers sell, in addition to the unit itself, basemaps for many places around the world together with a proprietary software program to upload the maps to the GPS. But with the vast and detailed coverage of OSM why not use this data as a background map in your GPS? Well, in fact you can, and it’s all doable with open source software.  There are tool chains available for both Windows and Linux.  The OSM Wiki has a detailed table of software for each OS.  I’m going to focus on Ubuntu Linux and two software packages: the java based mkgmap and GPS interface software QLandKarte GT .

On the OpenStreetMap website you can search for and download data using the “Export” tab, and choose “OpenStreetMap XML data” as the format. However you’re limited to rather small areas – usually less than 50 km. squared. If you know the longitude/latitude of the area you’re interested in, you can download larger regions directly from the Information Freeway website using wget. Here’s an example:

~/geodata/OSM$ wget -O bosnia.osm http://www.informationfreeway.org/api/0.6/map?bbox=17,42.5,19,44

The format after the “map?bbox=” is West,South,East,North in decimal degrees.

Now that you have your *.osm file, what’s the next step? I chose the mkgmap utility to process OSM files into a format that Garmin GPS can use. Ubuntu has this application in its repository so it was a quick

~/geodata/OSM$ sudo apt-get install mkgmap

to get it installed.  And it works fine with Ubuntu’s default Open-JDK. To use this tool you can either type all the options on the command line, or create a config file with the options entered, one per line. The configuration file looks like this:

~/geodata/OSM$ cat template.txt
mapname: 05042010
description: Bosnia&Herzegovina
country-name: Bosnia
country-abbr: BH
region: Sarajevo
gmapsupp
tdbfile
latin1

All the mkgmap options are well explained in the OSM Wiki. If you’re going to upload several OSM maps to the GPS  string them onto the command. And the tdbfile option will be required in the next step.  Now just run:

~/geodata/OSM$ mkgmap -c template.txt bosnia.osm

and you get a *.tdb file (along with the *.img files) which will be needed for upload to the GPS.

Next, to connect to my Garmin etrex Legend, I used the QLandKarte GT application.  While there’s an Ubuntu deb available for this software, it seems that it’s a minor version behind, and the latest version 0.18 has some fixes specifically for Garmin. So rather than use the binary from Ubuntu,  I downloaded both QLandKarte GT and their Garmin Drivers from their svn repository. For compiling, you’ll need cmake, qt4, and both proj4 and gdal,  all available from the regular Ubuntu and the UbuntuGIS repositories. You can compile following these instructions. I found one “gotcha” that needs to be fixed by hand. On Ubuntu 10.04 we already have gdal 1.7.2. But the cmake script that searches for GDAL only looks for versions up to 1.6.0. So I did a quick edit of the file:

<your download location>/QLandkarteGT/cmake/Modules/FindGDAL.cmake

and added the line:

gdal1.7.0

to the find_library(GDAL_LIBRARY)  section after line 47.

I didn’t need to change any other of the cmake defaults, and

make && sudo make install

ran successfully. Next I did cmake for the Garmin drivers, and also make && sudo make install, without any problems. It deserves mention that Ubuntu (and other distros) come with a kernel driver for garmin called garmin_gps, which apparently doesn’t work well. This kernel module needs to be removed and blacklisted. (Many distros by default have it blacklisted). So if you ask:

~/geodata/OSM$ lsmod | grep garmin
~/geodata/OSM$

If lsmod returns nothing, you’re OK. If, on the other hand, the module is loaded, then do

~/geodata/OSM$ sudo rmmod garmin_gps

and furthermore, edit your /etc/modprobe.d/blacklist.conf. Make sure it has a line

blacklist garmin_gps

We’re now ready to startup the QLandKarte application, and from the menu File->Load Map, choose the tdb file created earlier by mkgmap. You will then get a window asking which “basemap” to load. This will be the map name of 8 digits that you entered in the mkgmap configuration file earlier. The map will load (click Map->Center Map or F3 if it doesn’t appear). Next go into Setup->General->Device & Xfer to select your GPS make and model. Now from the Map menu choose Select Submap (or F5) and drag the mouse over the tiles you need.

On the bottom left should be the “Selected Maps” pane. Click on gmapsupp and then click the “Export” button. In the “Create gmapsupp” window first click the “browse” button on the right, and choose a destination directory on your computer, then click Export (creating this local copy makes the upload faster).

Screen shot of Garmin GPS with OSM layers as basemap

Finally, it’s time to fire up your GPS, and connect the USB cable to the computer. In the GPS Main menu, choose Setup->Interface. It should say USB (GARMIN format) and “Connected”. Now press the rocker button on the GPS, and you’ll see a blue screen with an image of a computer and GPS connected by a cable. You’re now in “Data Transfer” mode. Back in QLandKarte GT, in the Main window, click the Export button in the lower left corner again. But this time when you click the browse button go to the directory /media. You’ll see a directory named something like /media/6634-3422. That’s the Micro SD in your Garmin. You must put the gmapsupp.img into a subdirectory “garmin” on the Micro SD card.  If there’s no “garmin” directory, you should be able to open the Micro SD from Ubuntu’s “Places” and just create the dir. Then, from QLandKarte, use the Export button to place the gmapsupp.img file there.

Now close QLandKarte,  disconnect the GPS, and restart it. When the GPS comes back up you should have the new basemap available from the Setup->Map screen.

But suppose you need to manipulate attributes of a vector layer that’s already stored in a PostGIS spatial database. A typical example: I have a GIS layer of rain gauges, stored in PostGIS. I need to add a column to the layer’s attribute table for elevation, and then get the elevation for each gauge from a GRASS based DEM (Digital Elevation Model) raster.

There are three approaches. For this post I’ll call them:

1. Import-Export: The usual, long but straight-forward way
2. Duplicate table: Also somewhat long, but saves having two layers in PostGIS
3. Hop Skip and Jump: This suggestion I owe to Jonah Duckles. A bit harder to follow, but it avoids any duplication of layers or attribute tables.

The conventional Import-Export method involves importing the point vector from PostGIS into GRASS, adding the required attribute column and updating its values in GRASS, then exporting back out to PostGIS. But you end up with a second PostGIS vector layer, equivalent to the first with the addition of the required elevation column. Now you must use PostgreSQL tools to drop the original table and rename the new one to the original name.

Briefly…

(we’ll assume a PostGIS vector called gauges, and the GRASS elevation raster called dem)

# Set up the data source  details

export HOST=<Your PostGIS hostname>

export DB=<Your PostGIS database name>

export USER=<Your PostGIS administrator name>

export PASS=<Your PostGIS admin password>

export DSN="PG:host=${HOST} dbname=${DB} user=${USER} password=${PASS}"

v.in.ogr dsn=${DSN} layer=gauges out=gauges_tmp

# We now have a duplicate of the original PostGIS vector stored totally (geometry and attributes) in GRASS

# Add the new column for elevation values

v.db.addcol gauges_tmp col="elevation double precision"

# check attribute table

v.info -c gauges

# Now get the elevation values from the raster

v.what.rast vect=gauges_tmp rast=dem col=elevation

# check elevation values:

v.db.select gauges_tmp col=gauge_id,elevation

If everything looks good, export back to PostGIS:

v.out.ogr in=gauges_tmp dsn=$DSN type=point format=PostgreSQL

Now you should have in PostGIS two point layers, with the same features, and the gauges_tmp table should have the additional elevation column. What’s left is to drop (or rename) the original table and rename the new one to the original name. (Double check that the new gauges_tmp table contains all the features before continuing!)

psql -h ${HOST} -d ${DB} -U ${USER} -c "ALTER TABLE gauges RENAME TO gauges_orig"
psql -h ${HOST} -d ${DB} -U ${USER} -c "ALTER TABLE gauges_tmp RENAME TO gauges"

In PostGIS spatial tables also appear in the geometry_columns table so be sure to rename the gauges_tmp table to exactly the same name as the original.

Now on to the second Duplicate table method.
This way involves setting the database connection in GRASS to store attribute tables straight in PostgreSQL, then import as before the original gauges layer. This will create a GRASS vector in which the geometry is in GRASS (as always) but the attributes are kept right in PostgreSQL. Then run the v.what.rast command, and the new elevation values are updated right into the new PostGIS table. The final step is to copy those elevation values from this new attribute table to the original gauges layer. Here are the steps:

# First check current database connection and take note of the settings. We’ll need these parameters to revert to this database connection after the procedure.

db.connect -p >> connection.txt
db.connect driv=pg database="host=${HOST},dbname=${DB}"
db.login user=${USER} password=${PASS}

# and import the original gauges just like before

v.in.ogr dsn=${DSN} layer=gauges out=gauges_tmp

# This creates a GRASS vector, but the attribute table is now stored right in your PostGIS database. We’re ready to add the required elevation column and get values from the GRASS dem.

v.db.addcol gauges_tmp col="elevation double precision"
v.what.rast vect=gauges_tmp rast=dem col=elevation

All that remains now is to run two PostgreSQL commands: add an elevation column to the original table and update its values from the gauges_tmp table:

psql -h ${HOST} -d ${DB} -U ${USER} -c "ALTER TABLE gauges ADD COLUMN elevation float"
psql -h ${HOST} -d ${DB} -U ${USER} -c "UPDATE gauges SET elevation=gauges_tmp.elevation FROM gauges_tmp WHERE gauges_tmp.gauge_id=gauges.gauge_id"

The “WHERE” clause insures that the correct elevation values are entered into the permanent gauges layer based on the gauge_id column. Finish by removing the gauges_tmp layer (Its accompanying database attribute table will also be dropped).

g.remove vect=gauge_tmp

And don’t forget to reset the database connection back to the original values (We saved the parameters into the file connection.txt):

cat connection.txt
db.connect driv=... database=...

On to the Hop Skip and Jump method.

First I should note that this trick works only for point vectors. If you need to update attribute values for PostGIS based line or area features, this won’t work. Use one of the first two ways.

The process involves three routines that can be chained into a one-line command. First we take advantage of the GRASS v.external module to have access to externally stored vectors. This modules gives us a read-only representation in GRASS of an externally stored vector. Then, with the v.out.ascii command we get a list of the X-Y coordinates for all the points in this vector. These coordinates are piped as input to the r.what module which outputs the raster value at each of the input point locations. The output from this raster command is then fed into a psql statement to update the original database directly. No interim tables are needed at all. Here we go:

v.external dsn=${DSN} layer=gauges out=gauges_tmp

# Now the one-liner

v.out.ascii in=gauges_tmp fs=" " columns="gauge_id" | \
r.what in=dem fs=" " | \
sed '/\*/d' | \
awk '{print "UPDATE gauges SET elevation="$5 " WHERE gauge_id="$4}' >> update_gauges.sql

Here’s the “play by play”:

As I mentioned, v.out.ascii creates a list of X-Y coordinates, including the GRASS cat values, and also the gauge_id, added by the columns= parameter. This list is input to r.what which just adds at the end of each line the raster value at that X-Y location. But r.what will output an asterisk ‘*’ when no value is available (i.e for points outside the region of the raster). We don’t want these lines in the output, so the sed command deletes them.

The resulting list gets piped thru the awk command to create a list of PostgreSQL UPDATE commands, where the elevation values are taken from column 5 in the input, and the gauge_id is from column 4. (Columns 1,2 are the coords, column 3 the cat value. Thus column 4 is the gauge_id and column 5 the raster value).

# Check the results carefully for any anomolies…

less update_gauges.sql

# If all is OK, just run this file through psql:

psql -h ${HOST} -d ${DB} -U ${USER} -c "ALTER TABLE gauges ADD COLUMN elevation float"
psql -h ${HOST} -d ${DB} -U ${USER} -f update_gauges.sql

Quick and elegant. If you’re confident that the above string of commands creates correct sql UPDATE statements, you can even do away with the update_gauges.sql file and pipe the list straight into psql.

GM: I’m quite excited to be here.

MS: We have many questions for you, and much to discuss. Shall we begin with the map projection named after you, the Mercator projection? How did you come to realize that stretching the spaces between parallels would aid sea navigation?

GM: Actually it was the opposite. As exploration missions expanded and reached further abroad, we were getting measurements back from navigators that were inconsistent. The more we investigated, we began to notice that navigators would mark the positions of coastlines that they reached on their captain’s logs after sailing for weeks at a constant compass direction. But a rhumb line on a sphere is not a straight line on a flat map. So we corrected that with our new representation of the earth.

MS: Did you understand at the time the severe distortion in area and size that your map projection introduced?

GM: I think you exaggerate. After all, distortion appears only near the poles in the far north and south. We had very little interest in these areas. What’s more, when sailing north or south, near to the meridian of departure, sea navigators could find there location quite well with the mariner’s quadrant. Only those explorers who ventured far west or east encountered difficulty in determining their location since they had no reliable chronometer.

MS: Why was that so important?

GM: Well, I’m surprised you ask.  As any beginning student of geography understands, a navigator can find his latitude very accurately by siting the angle of Polaris. However the only way he can determine his longitude is by comparing the angle of the sun or stars in the sky at his current position with the angle at some know location at exactly the same time of day. So he must know the hour precisely. And, of course, we have no instrument which could keep time accurately at sea, with the swaying of the vessel, and corrosive effects of salt spray. Perhaps some new, advanced technologies have led to a reliable, sea-worthy chronograph?

MS: Yes, well the first real chronometer for navigation was build by an English clock maker, John Harrison, 100 years after your death. Today we have clocks and watches that are water-proof, and stay very accurate for years. But navigation as well as map making now uses GPS.

GM: GP what?

MS: There are 24 satellites in orbit around the earth at about 11,000 km altitude. They circle the earth twice every day.

GM: Mijn God! You have man made instruments in the heavens! And they imprisoned  me as a heretic for drawing different maps… So what do these “satellites” do?

MS: They constantly send a signal which contains the satellite’s identity and the precise time.  Anyone using a simple hand-held GPS can receive those signals, and the instrument can determine your location to within 10 meters.

GM: This all sounds like a fairy tale. How could such a small instrument possibly find locations so accurately?

MS: As I said, the instrument picks up signals from several of the satellites, and it compares the time differences from the different satellites.  From the gaps in time it calculates the distance to each of the satellites, and thus by doing a kind of triangulation of these distances, it can zero in on the current location.

GM: I shall not believe any of this until I see it for myself.  In any case, map makers today will still need my method to properly spread the terrestrial globe unto a planar map.

MS: In fact we have thousands of map projections in use today. Each country has it’s own local grid. Some have several regional grids. All are designed to ensure that surveys done in different locations are all synchronized and match. So, you see that your dream of creating a world map by combining many local, independent maps has, in a round about way, been fulfilled.

GM: So many different methods of projection. That must be a nightmare to transform surveys from one projection to another.

MS: That’s all done by computers.

GM: You must employ many of these computers, and they all must be highly trained in mathematics, earning respectable salaries.

MS: No, computers are electronic instruments that perform mathematical calculations very rapidly. These computers are programmed with instructions we call “software” to carry out many, many complex calculations. In fact, today cartographers use computers to design and produce maps. All the geographical data is kept inside the computer’s memory, and the cartographer can choose which data he wants, and how to display it. Many versions of a map, at many different scales can be produced and printed in a short time.

GM: So you have instruments orbiting in the heavens, instruments for finding your location, instruments for doing mathematics and creating maps. This is all very fatiguing.  I must take a rest now.

So, suppose we want to make a map displaying both Lon-Lat geographic coordinates, and a local CRS. Begin by loading a layer covering the region of interest and setting the project CRS to EPSG:4326 (Lon-Lat, WGS84) with on-the-flye projection enabled. Now use the “Vector->Research Tools->Vector Grid” option to create a grid with the spacing and extents required. Once the grid shapefile is created, open it’s *.dbf file in Openoffice.org Calc and, following the details in the above post, add columns for X_OFFSET, Y_OFFSET, ANGLE, and INTERVAL. The offsets you will enter refer to the current CRS (Lon-Lat WGS84 in this case). Once you have that setup, display this new grid shapefile with labels. You’ll use the label options (as explained in the same post mentioned above) to set the angle and X and Y coordinates to get the labels placed around the edges of your map.

Now change the project’s CRS to whatever the second grid should be, and maintain on-the-fly projection enabled. You’ll notice that the first Lon-Lat grid will now be slightly “distorted” due to the reprojection.  Open the Print Composer, and add a Map element. In the Map Item configuration window, go to the Grid tab, and check “Show Grid?”. Set up the intervals, and fonts, etc for the grid. This grid will represent the project’s projected coordinate system.

Here’s an example. The data were obtained from the Natural Earth data website, and reprojected within QGIS to Jordan Transverse Mercator – EPSG:3066 (the blue lines). The red coordinate grid represents Lon-Lat geographic coordinates:

Map of Middle east, Jordan Transverse Mercator projection

To follow up my last post on isohyetal lines, I want to check out the new Contours plugin authored by Lionel Roubeyrie and compare the results to contour lines created with GRASS. With no futhur ado here’s the result:

Two methods for creating contour lines

The green lines are my original isohyetal lines (made with GRASS). The blue lines are based on the same rainfall data  with the Contour plugin, and choosing 12 levels of contours. The plugin doesn’t have an option to set interval, rather the number of levels. Thus the values will not be identical to the original run. In the center area the two sets of lines pretty much follow each other. What is obvious to me is the “confusion” around the edges. While the GRASS algorithm does fairly good extrapolation outside the area covered by data, the Contour plugin seems to just “cut off”. So the next test should be with a higher density of points and restricted to the region actually covered by data.

]]> http://www.surfaces.co.il/?feed=rss2&p=595 5 http://creativecommons.org/licenses/by-nc-sa/3.0/ http://www.surfaces.co.il/?p=578 http://www.surfaces.co.il/?p=578#comments Wed, 10 Feb 2010 19:01:24 +0000 Micha Silver http://www.surfaces.co.il/?p=578 We had a pretty extraordinary rain storm in our region some weeks ago. Accumulated rainfall over a 24 hr period was between 20-100 mm in a region where the total annual precipitation is about 50 mm.! I got rain gauge data for the event and made isohyetal lines using only tools available in Quantum GIS.

First I imported the data table including XY locations and the precipitation with the Delimited Text plugin. Then, in order to limit the analysis more or less to the area covered by the rain gauges, I made a “convex hull”, the minimum polygon enclosing all points, and buffered that polygon by 10 km. (allowing that the interpolation algorithm will give approximate values outside the area covered by the gauges).  Both of these operations are available in the Vector->Geoprocessing menu.

Now I fired up the GRASS plugin to do the interpolation. Using the v.in.ogr.qgis module, I loaded both the rain gauge point vector and the buffer polygon vector into a suitable GRASS Location/Mapset. First I converted the buffer polygon to a raster so that I could use it as a mask with v.to.rast. After adding the two GRASS layers to the map – the rain gauges, and the mask raster,  I pulled up the r.mask module to force the next action to be limited to the buffer region. Then I ran v.surf.rst to produce an interpolated rainfall grid. I chose, of course, the precipitation column as the attribute field for doing the interpolation. The new precipitation grid was created in a few moments, and I closed the GRASS toolbox.

Now I activated the new GDALTools Raster plugin. Among the tools there is “Contours”. I ran this tool, choosing the GRASS precipitation raster as input. I left the default levels value at 10, and chose an output directory where the contours shapefile will be saved. I also checked “Attribute Name” and typed in “Precip”. The contours were created and here’s my resulting map:

Precipitation map - january 2010

GRASS

Now we’re ready (almost) for compiling GRASS. First get more pre-requisites:

~/Download$ sudo yum install tcl tcl-devel tk-devel fftw2 fftw2-devel libXmu libXmu-devel

The Swig conundrum:  SWIG is an interface that “dresses” C/C++ programs for use in scripting languages like python. It seems that the installed version in F12, 1.3.40, has some change that makes it incompatible with certain parts of the GRASS python interface. The recommendation at this point is to roll back to an earlier version of SWIG (Note: if you skip this step, GRASS will be fully functional in both the CLI and the tcl/tk interface but some parts of the wxpython interface won’t work). You’ll have to remove the rpm and download the older 1.3.36 version. To successfully do the compile you need the perl devel package (and others that we’ve already added). So here are the steps:

~/Download$ sudo yum install perl-devel
~/Download$ sudo yum remove swig
~/Download$ tar xvzf swig-1.3.36.tar.gz
~/Download$ cd swig-1.3.36
~/Download/swig-1.3.36$ ./configure
~/Download/swig-1.3.36$ make && make check
~/Download/swig-1.3.36$ sudo make install

On to GRASS itself. Download the  GRASS source code for version 6.4, and as always, unzip the source and go into that directory.

Here’s my configure line:

~/Download/grass-6.4.0RC5$ ./configure --prefix=/usr/local/grass-6.4  \
      --with-proj-share=/usr/share/proj  \
      --with-postgres  \
      --enable-largefile  \
      --with-sqlite  \
      --with-freetype --with-freetype-includes=/usr/include/freetype2  \
      --with-python --with-wxwidgets  \
      --with-cxx

The freetype include headers are in an unusual place: /usr/include/freetype2/freetype. GRASS apparently always looks in a subdir called ‘include’.  So in order to get GRASS to recognize the freetype headers I had to do:

~/Download/grass-6.4.0RC5$ cd /usr/include/freetype2
/usr/include/freetype2$ sudo ln -s freetype include

then the above compile line worked. Also don’t forget to add –with-cxx, (otherwise you won’t have the python digitizer)

Now run make and make install

/usr/include/freetype2$ cd ~/Download/grass-6.4.0RC5
~/Download/grass-6.4.0RC5$ make
~/Download/grass-6.4.0RC5$ sudo make install

Next we’ll be installing Quantum GIS. The path to the grass lib directory must be added to a conf file in ld.so.conf.d (for the same reason explained in the earlier post in the ECW section). In my case I create a file  /etc/ld.so.conf.d/grass-6.4 which contains one line:

/usr/local/grass-6.4/grass-6.4.0RC5/lib

then run: $sudo ldconfig

At this point we have one more step before moving on to Qunatum GIS. We compiled GDAL earlier without GRASS support.  That’s fine if you don’t need access to GRASS rasters when using the gdal tools. If you do want to manipulate GRASS rasters using GDAL tools, or more importantly, in order to access GRASS raster from QGIS (which depends on GDAL) then you want the gdal-grass plugin. It’s quick and simple. Download the lastest version 1.4.3 and after opening the tar archive do the usual:

~/Download/gdal-grass-1.4.3$ ./configure --with-grass=/usr/local/grass-6.4/grass-6.4.0RC5
~/Download/gdal-grass-1.4.3$ make && sudo make install

QGIS

Prepare for Quantum GIS with:

~/Download$ sudo yum install cmake gsl gsl-devel geos-devel expat expat-devel

Download the QGIS source code for version 1.4, unzip it and drop into that directory.

Following more or less Tim Sutton’s blog post or my earlier instructions regarding Fedora 10, invoke these commands:

~/Download$ cd qgis-1.4.0/
~/Download$ mkdir build
~/Download/qgis-1.4.0$ cd build
~/Download/qgis-1.4.0/build$ cmake -L ..

(Be sure to add the ‘..’ at the end of the cmake command).  Now check thru the list of cmake parameters to make sure all the dependencies were found.  You tune these paths or variables as necessary with the -D option to cmake. In my case I need to indicate where GRASS is located (note the –prefix option in my GRASS compilation previously), and I want to set the final install location to a specific directory under /usr/local for this version of QGIS, so:

~/Download/qgis-1.4.0/build$ cmake -L -DGRASS_PREFIX:PATH=/usr/local/grass-6.4/grass-6.4.0RC5 -DCMAKE_INSTALL_PREFIX:PATH=/usr/local/qgis-1.4.0 ..

With that set up we’re ready for:

~/Download/qgis-1.4.0/build$ make
~/Download/qgis-1.4.0/build$ sudo make install

R statistics

No GIS workstation is complete without  R and its spatial packages.

~/Download$ sudo yum install R-core rpy R-devel xml2 libxml2 libxml2-devel

Fire up R and run:

> install.packages("sp",dependencies=TRUE)
> install.packages("spgrass6",dependencies=TRUE)

The R install.packages() function wants to put R libraries into /usr/lib/R/library and the html manuals into /usr/share/doc/R-x.x.x by default. Since these locations are writable only by root, either give yourself necessary permissions on those directories, or (preferably) in the install.packages() functions add the lib=”…” parameter to create an R library in your home dir.

Whew! Congratulations, you’re now the proud owner of an advanced Linux based GIS workstation.

I won’t go into the details of wading thru the search and order procedure. The site has it’s own  tutorial. Once you complete the search, and get a list of “granules” of data that match your search criteria, you tick which you want, and click the “Add to cart” button. Then enter some required order details, and click “Submit order”. There are buttons to save query and personal order details to your computer. I suggest to use this feature: makes the next visit smoother.

Some minutes later you’ll receive three emails. The second, with a subject like: “LPDAAC ECS Order Notification” contains the ftp link to download the data you requested. I copy the link and use wget to pull down the zip file.

~/geodata/MODIS_AQUA$ wget -c ftp://e4ftl01u.ecs.nasa.gov/PullDir/0301944518RAPMac.zip

After unzipping the file, l have several *.hdf and *.xml files. Now I can run the hdp utility to quickly check that the file is OK.

~/geodata/MODIS_AQUA$ /usr/local/hdf4.2/bin/hdp dumpsds -h MYD13A3.A2003032.h20v04.005.2008056162517.hdf

which should output a dump of all the headers and data in the file.  Now we’re ready for bigger and better things.

GDAL

Compile and install GDAL. The source is available here for version 1.6.3. Unzip and go into the source directory.

~/Download$ cd gdal-1.6.3

first some additional pre-requisites

~/Download/gdal-1.6.3$ sudo yum install sqlite sqlite-devel
~/Download/gdal-1.6.3$ sudo yum install postgresql postgresql-devel postgresql-server postgis pgadmin3

It seems that the external libgeotiff include header files are in their own subdirectory.  So here’s the configure incantation I used:

~/Download/gdal-1.6.3$ CPPFLAGS="-I/usr/include/libgeotiff" ./configure --with-geos --with-sqlite3 --with-ecw=/usr/local/ --with-hdf4=/usr/local/hdf4.2 --with-python
~/Download/gdal-1.6.3$ make && sudo make install
~/Download/gdal-1.6.3$ sudo ldconfig

Now check that the formats we compiled in are indeed available with:

~/Download/gdal-1.6.3$ gdalinfo --formats | egrep 'HDF|ECW'
  HDF4 (ro): Hierarchical Data Format Release 4
  HDF4Image (rw+): HDF4 Dataset
  ECW (rw): ERMapper Compressed Wavelets
  JP2ECW (rw+): ERMapper JPEG2000

Good. We can now put gdal to work. First, I run gdalinfo on one of the MODIS hdf files.

~/geodata/MODIS_AQUA$ gdalinfo MYD13A3.A2003032.h20v04.005.2008056162517.hdf

This gives me several lines that start with SUBDATASET_x_NAME (‘x is the number of the dataset).  Let’s run gdal again and filter the results with egrep:

~/geodata/MODIS_AQUA$ gdalinfo MYD13A3.A2003032.h20v04.005.2008056162517.hdf | egrep "SUBDATASET_[0-9]_NAME"
 SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MYD13A3.A2003032.h20v04.005.2008056162517.hdf":MOD_Grid_monthly_1km_VI:1 km monthly NDVI
 SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"MYD13A3.A2003032.h20v04.005.2008056162517.hdf":MOD_Grid_monthly_1km_VI:1 km monthly EVI
 SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MYD13A3.A2003032.h20v04.005.2008056162517.hdf":MOD_Grid_monthly_1km_VI:1 km monthly VI Quality
 ....

Looking thru the results I see that SUBDATASET_2 is the EVI data. We can now run gdalinfo on that single SUBDATASET by feeding it everything on the line after the “…_NAME=” like so:

~/geodata/MODIS_AQUA/mesopotamia$ gdalinfo 'HDF4_EOS:EOS_GRID:"MYD13A3.A2003032.h20v04.005.2008056162517.hdf":MOD_Grid_monthly_1km_VI:1 km monthly EVI'

Note that the whole dataset string is ‘quoted’. The output now shows all the georeference information that gdal finds in that particular  dataset. Next step is to reach into the gdal toolbox again, and pull out gdalwarp. That will allow us to both export the data to a GeoTiff and at the same time to transform it to a CRS that we choose. The input to gdalwarp must be the full dataset name as in the last run of gdalinfo above.

~/geodata/MODIS_AQUA$ gdalwarp -of GTiff -t_srs EPSG:4326 -dstnodata "-3000" -dstalpha 'HDF4_EOS:EOS_GRID:"MYD13A3.A2003032.h20v04.005.2008056162517.hdf":MOD_Grid_monthly_1km_VI:1 km monthly EVI' EVI_h20v04.tif

Note that I add the -dstnodata and dstalpha options to gdalwarp to set the “nodata” value in the EVI to NULL. Why “-3000″ ?  The EVI dataset uses -3000 as it’s nodata value. Furthermore, a multiplier of 10,000 has been applied so that all values in the HDF are integer. More on that in the next post.

If you’ve downloaded several HDF files, and converted them to GeoTiff as above, then you probably want to merge them together into one. GDAL offers the gdal_merge.py utility for that:

~/geodata/MODIS_AQUA$ gdal_merge.py -o EVI_mesop.tif EVI_*.tif

I downloaded vegetation index files for the Tigris-Euphrates rivers drainage basin, thus I chose to name the merged tif “Mesopotamia”.

HDF Java tools

With that under our belt, it’s time for me to come clean: there’s an easier way to deal with HDF files. The hdfgroup has some java based utilities that make life simpler: First I’ll mention HDFView to select and view datasets within an HDF file. And, more importantly, NASA offers HEG, the HDF-EOS to GeoTiff conversion tool. This tool does what we did above with gdal: converts to GeoTiff, transforms to a coordinate reference system, and as an extra perk, it merges several tiles at the same time. The installation of both tools is a java binary. So once you’ve downloaded them:

~/Download$ sudo sh ./hdfview_install_linux32.bin

The java installer is straightforward. Then run the tool  (if you installled in /usr/local) by typing: /usr/local/hdfview/hdfview.  In the java interface you can open an HDF file and drill down to a specific dataset. Then you right-click on the dataset, and choose “Open As” and select Image to display an image of the data.

The HEG tool is somewhat tricky. First, it requires full write permission on the install directory. So if you are sticking with /usr/local, you must give yourself write permissions there. (Alternatively, just install in your home directory). Second, the tool does not recognize any $PATH variables. So it must be run from the install directory.

micha@hayun-13:~/Download$ tar xvzf hegLNXv2.10.tar.gz
install
heg.tar
micha@hayun-13:~/Download$ sudo ./install
micha@hayun-13:~/Download$ sudo chown -R micha /usr/local/heg

As part of the installation, you enter the location of the java binary, /usr/bin/java (Fedora’s “/etc/alternatives” infrastructure takes care of that), and the final installation location, /usr/local/heg. Once the installation is complete, run the program by:

micha@hayun-13:~/Download$ cd /usr/local/heg/bin/
micha@hayun-13:/usr/local/heg/bin$ ./HEG

The java interface has it’s own Help. Select a file from the File->Open menu, and a dataset in the left pane. Then the output file name and projection parameters in the center pane. Click “Accept”, and “Run” in the right hand pane. If, before clicking Run, you open other files, select a dataset, and click “Accept”, then the tiles are merged into one big GeoTiff .

The above map was put together using the techniques above. The EVI data is from the MODIS AQUA dataset from May 2009, Monthly average vegetation Index. The vector layers come from the public domain dataset Natural Earth. I used a few GRASS operations to merge several tiles and get the coloring right, but that’s for the next blog… So we’ve concluded installing the foundations for GIS software. Tune in next time for compiling QGIS and GRASS.