Category Archives: GIS

Fourth release candidate of GRASS 6.4.0 now available

A fourth release candidate of GRASS 6.4.0 is now available.

Source code:

To get the RC4 source code from SVN:
svn checkout

An announcement has been drafted at
(all RC news will be merged into the final announcement later)

Key improvements of the GRASS 6.4.0 release include enhanced portability for MS-Windows (native support), hundreds of fixes, the new wxPython based portable graphical interface and much new functionality.

Release candidate management at

We hope to get out binary packages over the next days.

Top ten myths for Open Source in Geospatial

An interesting email thread is ongoing about Top ten myths which try to harm the reputation of Open Source GIS efforts. Michael P. Gerlek edits a Wiki page collecting the top ~10 myths and misperceptions:

Especially striking the Eduardo Kanegae‘s comment:

In 2008 I worked on a project all based on ESRI 9.2 family. At that
point I didn´t know much of ESRI products and had only worked with
foss products. Now I feel more confortable to give an opinion for

# myth monopoly : every only remember the (supposed) 30% esri market
share. Remember there´s also very nice commercial products like Safe
FME, CadCorp, ManiFold, PCI, ERDAS, ENVI and others ( and sometimes
most of these are embed on ESRI packs - eg.: raster support on AG
family )

# sustainability : while every major release of ESRI will force you to
re-develop your customizations, FOSS products keep release more
compatible. Example: a MapServer 3.x developer will use the same
principles and concepts on MapServer 5.x version. But, an ArcIMS
developer had to change its base when upgrading to ArcGIS Server 9.1,
recode for adapting to ArcGIS Server 9.2 API´s and now all this
concepts will change again with 9.3 version.

# maintenance : foss product will run more closer to open standards (
eg.: OGC´s ). So, you change foss parts without re-coding your entire
solution. The cost of training a new human resource on
insert/update/delete geo-feature using ArcObjects/ArcSDE is so much
higher when compared to OGC-SFS, per example.

# support : while on FOSS communities you can have a reply on minutes,
on 'esri forum' you can your topic open for months ( , , ) and NEVER
get a solution. In my sample, we discover a bug on ArcSDE/ArcMap
9.2sp4 but this will certainly NEVER be fixed because sp6 didn´t fixed
and ESRI will probably only look for 9.3 developments for now and on.
Because non-US customers CAN´T contact ESRI directly, we can only keep
suffering with poor local support.

best regards,
Eduardo Kanegae

Not much to add at this point…

OGR vector data tips and tricks

Get vector map extent

You can easily grep the map extent of a vector map (bounding box):

ogrinfo /cdrom/ITALY_GC.SHP ITALY_GC | grep Extent
Extent: (9.299460, 43.787362) – (13.911350, 47.070188)

Merge of two SHAPE files

Merge of two SHAPE files ‘file1.shp’ and ‘file2.shp’ into a new file ‘file_merged.shp’ is performed like this:

ogr2ogr file_merged.shp file1.shp
ogr2ogr -update -append file_merged.shp file2.shp -nln file_merged file2

The second command is opening file_merged.shp in update mode, and trying to find existing layers and append the features being copied. The -nln option sets the name of the layer to be copied to.

Vector map reprojection

We reproject from the source projection (as defined in .prj file) to WGS84/LL:

ogr2ogr vmap0rd_ll.shp -t_srs “EPSG:4326” vmap0rd.shp

If the .prj file is missing, you can use the ‘’ utility to create it if you know the EPSG code: -wkt 4326 > cities.prj

Reproject to current GRASS location projection:

ogr2ogr -t_srs “`g.proj -wf`” polbnda_italy_GB_ovest.shp polbnda_italy_LL.shp

Cut out a piece of a vector map
Use spatial query extents: -spat xmin ymin xmax ymax (W S E N)

ogr2ogr ARC_BZ.shp -spat 10 45 13 47 ARC.shp

Get VMAP0 metadata info:

ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr | grep bnd
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘roadl@trans(*)_line’

MAP0: Extract spatial subregion, reproject from NAD83 to WGS84

# coordinate order: W S E N
ogr2ogr -spat 19.95035 -26.94755 29.42989 -17.72624 -t_srs ‘EPSG:4326’ \
polbnda_botswana.shp gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr \


Sample ‘where’ statements (use -sql for PostgreSQL driver):

# -where ‘fac_id in (195,196)’
# -where ‘fac_id = 195’
ogrinfo -ro -where ‘fac_id in (195,196)’ \
gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’

VMAP0 examples

Find out the Countries VMAP0 coding:

ogdi_info -u gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr \
-l ‘polbnda@bnd(*)_area’ -f area | grep Botswana

or read the VMAP0 Military specs, page 75

Extract Botswana, reproject on the fly from NAD83 to WGS84, store as SHAPE:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘BC'” polbnda_botswana.shp \
gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’

Extract Germany, reproject on the fly from NAD83 to WGS84, store as SHAPE:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘GM'” polbnda_germany.shp
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘polbnda@bnd(*)_area’
ogrinfo -summary polbnda_germany.shp polbnda_germany | grep Extent
# Extent: (5.865639, 47.275776) – (15.039889, 55.055637)
# W S E N

VMAP0 Contour lines for Germany:

ogr2ogr -t_srs “EPSG:4326” -spat 5.865639 47.275776 15.039889 55.055637 \
contour_lines.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘contourl@elev(*)_line’

VMAP0 elevation spots (points) for Germany:

ogr2ogr -t_srs “EPSG:4326” -spat 5.865639 47.275776 15.039889 55.055637 \
elevation_spots.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘elevp@elev(*)_point’

VMAP0 lakes of Trentino province in Italy:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘IT'” \
-spat 10.340029 45.261888 10.98727 45.98993 \
lakes_italy.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘inwatera@hydro(*)_area’

Connect OGR and PostgreSQL/PostGIS

ogrinfo PG:’ user=postgres dbname=ogc_simple’
ogr2ogr out.shape PG:’ user=postgres dbname=ogc_simple’ lake_geom


Convert GRASS 6 vector map to SHAPE (needs GDAL-OGR-GRASS plugin):

# -nln is “new layer name” for the result:
ogr2ogr archsites.shp grassdata/spearfish60/PERMANENT/vector/archsites/head 1 \
-nln archsites

Using WKT files with ogr2ogr

The definition is in ESRI WKT format. If you save it to a text file called out.wkt you can do the following in a translation to reproject input latlong points to this coordinate system:

ogr2ogr -s_srs WGS84 -t_srs ESRI::out.wkt out_dir indatasource

Most comand line options for GDAL/OGR tools that accept a coordinate system will allow you to give the name of a file containing WKT. And if you prefix the filename with ESRI:: the library will interprete the WKT as being ESRI WKT and convert to “standard” format accordingly. The -s_srs switch is assigning a source coordinate system to your input data (in case it didn’t have this properly defined already), and the -t_srs is defining a target coordinate system to reproject to.

TIGER files in OGR

# linear features:
ogr2ogr tiger_lines.shp tgr46081.rt1 CompleteChain

# area features:
export PYTHONPATH=/usr/local/lib/python2.5/site-packages tgr46081.rt1 tiger_area.shp

OGR CSV driver: easily indicate column types

You can now write a little csv help file to indicate the columns types to OGR. It works as follows. Suppose you have a foobar.csv file that looks like this:


Now write a foobar.csvt file like this one:


The driver will then use the types you specified for the csv columns. The types recognized are Integer, Real and String, DateTime, and Date.

Convert KML to CSV (WKT)

First find layers:

ogrinfo -so myfile.kml

Then convert KML to CSV:

ogr2ogr -f CSV out.csv myfile.kml -sql “select *,OGR_GEOM_WKT from myfilelayer”
cat out.csv

Or use the cool online converter:

GDAL raster data tips and tricks

Reading GRASS data through GDAL/OGR support

Example 1: We write out a GRASS raster map to GeoTIFF — this format
includes the coordinates within the file’s metadata:

gdal_translate -of Gtiff /usr/local/share/grassdata/spearfish/PERMANENT/cellhd/soils soilmap.tif

ogr2ogr roadsmap.shp /usr/local/share/grassdata/spearfish/PERMANENT/vector/roads/head

Fast image display with tiling
If you want fast access you might want to try converting e.g. a BIL files to a tiled TIFF, and build overviews. You can build overviews for BIL too, but it can’t be directly tiled:

gdal_translate source_bil global30.tif -co “TILED=YES” -co “TFW=YES” -co “PROFILE=BASELINE”
gdaladdo global30.tif 2 4 8 16

GDAL performance problem?
GDAL_CACHEMAX is normally a number of megabytes (default is 10 or so). So something like:
gdal_translate -of GTIFF -co TILED=YES –config GDAL_CACHEMAX 120 madison_1f_01.jpg madison_1f_01.tif
would use a 120MB cache.

GDAL and 1 bit maps
With a trick you can get those: -co NBITS=1 -o dst.tif src.tif

Generate 8 bit maps for Mapserver
gdal_translate -scale in.tif out.tif
Note: As of MapServer 4.4 support has been added for classifying non-8bit raster inputs

Greyscale conversion
A “proper” conversion would involve a colorspace transformation on the RGB image into IHS or something like that, and then taking the intensity. GRASS can do things like that.

Generate an OGC WKT (SRS)
In WKT the ellipsoid is described by two parameters: the semi-major axis and the inverse flattening. For a sphere the flattening is 0 and so the inverse flattening is infinity.

# in the GDAL source code:
cd apps
make testepsg

./testepsg ‘+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m’
Validate Succeeds.
WKT[+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m] =
GEOGCS[“GRS 1980(IUGG, 1980)”,

Simplified WKT[+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m] =
GEOGCS[“GRS 1980(IUGG, 1980)”,

Extracting spatial subset (subregion)
gdal_translate -of GTiff -projwin 636861 5152686 745617 5054047.5 \
p192r28_5t19920809_nn1.tif test1_utm.tif

Fixing broken projection/datum info for raster data
gdal_translate -of HFA -a_srs epsg:32735 /cdrom/173072lsat.img \

# or, using a WKT file
gdal_translate -of HFA -a_srs file.prj /cdrom/173072lsat.img \

Merge various import maps, re-project on the fly and extract spatial subset according to current GRASS region
eval `g.region -g`
gdalwarp -te $w $s $e $n *.TIF \

Export to (limited) TIFF readers such as ArcView, or ImageMagick
Many tools have trouble reading multi-band TIFFs with “band interleaving”, the GDAL output default. Best is to use the INTERLEAVE=PIXEL creation option. Just add to the gdal_translate command line:

Inserting metadata (metadata tags)
gdal_translate -outsize 37.5% 37.5% \
in.tif out.tif

Raster map reprojection (warping)
gdalwarp -t_srs ‘+init=epsg:26591 +towgs84=-225,-65,9’ test1.tif \

Raster map reprojection (warping) maintaining NULL values (sea etc):

gdalwarp -r bilinear -tr 1000 1000 \
-srcnodata “-32768” -dstnodata “-32768” \
-wo “INIT_DEST=-32768” \

-t_srs epsg:32632 italy_LL.tif italy_UTM32.tif

Reprojecting external map to current GRASS location externally
gdalwarp -t_srs “`g.proj -wf`” aster.tif aster_tmerc.tif

Cut out region of interest with gdalwarp (in target coords)
Add to command line (insert values instead of letters of course:
#damn order, differs from -projwin!!
-te W S E N

Merging many small adjacent DEMs into one big map (A)
This needs GDAL compiled with Python and numpy installed:
# if not installed in standard site-packages directory
export PYTHONPATH=/usr/local/lib/python2.5/site-packages -v -o spearfishdem.tif -n “-32768” d*.tif

Merging many small adjacent DEMs into one big map (B)
Even easier, just use gdalwarp:
gdalwarp C_1mX1m/dtm*.tif big.tif
Or just a few tiles:
gdalwarp C_1mX1m/dtm0010[4-5]* big_selection.tif

Merge various map/bands into a RGB composite -of HFA -separate band1.img band2.img band3.img -o out.img

GDAL gdalwarp interpolation comments
Which method -rn, rb, -rc or -rcs should one use for DEM and which for data like e.g. Landsat TM reprojecting?

-tps: Enable use of thin plate spline transformer based on available GCPs.
-rn: Use nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
-rb: Use bilinear resampling.
-rc: Use cubic resampling.
-rcs: Use cubic spline resampling (slowest algorithm).

FrankW suggests:
I would suggest -rb for DEMs, and one of the cubic kernels for landsat data. Of course, there are various factors that you should take into account. Using -rb (bilinear) for the DEM will perform local averaging of the nearby pixel values in the source. This give reasonable results without introducing any risky “overshoot” effects you might see with cubic that could be disturbing for analysis or visualization in a DEM. The cubic should in theory do better at preserving edges and general visual crispness than using bilinar or nearest neighbour. However, if you are wanting to do analysis with the landsat (such as multispectral classification) I would suggest just using -rn (nearest neighbour) so as to avoid causing odd effects to the spectral values.
Nobody can’t tell you what method should be used in your case. Generally speaking, in the case of upsampling spline and cubic interpolators are more suitable (-rcs and -rc). In the case of downsampling and the same resolution it is completely up to you what method looks better. Just try them all and select the one which is most appropriate for you.

Geocoding with ‘gdal_translate’
FrankW suggests:
As far as I know there is not on-screen method for doing this, but it certainly isn’t too difficult with a little bit of semi-manual work. Open two OpenEV views, one with the unreferenced image, one with the geo-reference base you want to use. Move your cursor over the non-referenced one (let’s call it image1), record (read: write down!) the pixel x,y values. Then look at the same location in image2. Write down the geocoordinate for the pixel. You should have four numbers for each location you want to pin the image to. And so on and so on. Then use gdal_translate to translate image1.tif to image1_georefd.tif but adding the -GCP parameter for each set of coordinates. Like so…

gdal_translate -gcp 1 1 500000 5000000 \
-gcp 200 400 550000 5250000 image1.tif \

gdalinfo pg-PR1B0000-2002031402_100_001

To select a channel and warp to UTM (or whatever is inside):

gdalwarp HDF4_SDS:ASTER_L1B:”pg-PR1B0000-2002031402_100_001″:2 aster_2.tif
gdalinfo aster_2.tif

Impressive GRASS GIS Web site statistics

Just to update you on the GRASS Web statistics development, here the statistics (remember, we have MANY mirror sites):

Month Unique Number Pages Hits Bandwidth
visitors of
Jan 2008 39223 74088 291166 715946 101.23 GB
Feb 2008 38984 74043 218314 623770 107.09 GB
Mar 2008 40674 73389 223666 621816 107.04 GB
Apr 2008 5490 15702 135134 403726 220.87 GB
May 2008 20613 104556 912263 2242942 1442.31 GB

(this includes of course search engine traffic)

It appears that many visitors came back in May who downloaded the long awaited GRASS 6.3.0 release from 23 Apr 2008.

Some outstanding hits for May (views, only
10095 /grass63/binary/mswindows/native/
3271 /grass63/binary/mswindows/native/WinGRASS-6.3.0-Setup.exe

This points out of obvious need for a portable, in this case also MS-Windows compliant GIS which GRASS 6.3.0 now is! Fetch native winGRASS with installer or GRASS for MacOSX or GRASS for Linux or …

d.vect: support for z height (geometry) colors added

After getting mad with Lidar points colorizing which till now
required a DB table with GRASSRGB attributes, I have
modified d.vect to support colors directly from z height (geometry).
Works for 3D points, lines (eg, 3D contours) and 3D polygons
(eg delaunay triangles):

# Spearfish:
g.region rast=elevation.10m
r.random elevation.10m n=5000 vector=random3d -d
d.mon x0
# display as black points
d.vect random3d
# display 3D points colorized according to z height
d.vect -z random3d zcol=gyr

# 3D contour lines
r.contour elevation.10m out=contour20m step=20
d.vect -z contour20m zcol=gyr

# generate 3D triangles
v.delaunay random3d out=random3d_del
# display 3D polygons colorized according to z height
d.vect -z random3d_del type=area zcol=gyr


Proposal: OSGeo Cartographic Library

GRASS, QGIS and others are in the need of their own map printing tools for high quality output but these projects should not start from scratch. There is a wealth of underlying code already in Mapserver, Mapguide etc which could be re-used in the terms of their respective licenses and certainly of programming language compatibility. A joint OSGeo Cartographic Library comes to mind.

See also the ongoing related discussion!

GRASS GIS releases version 6.3.0

GRASS GIS releases version 6.3.0
23 April 2008

GRASS 6.3.0 is a “technology preview” release, the first beta on the path to GRASS 6.4-stable, and also marks the start of work on GRASS 7. As such GRASS 6.3.0 is not intended to be a stable release with ongoing support, but after five months of quality-assurance review users can be confident to use this version for their day to day work, indeed due to the open development model many already do.

This release brings hundreds of new module features, supported data formats, and language translations, as well as a number of exciting enhancements to the GIS. A prototype of the new wxPython user interface is debuted, and for the the first since its inception with a port from the VAX 11/780 in 1983, GRASS will run on a non-UNIX based platform: MS-Windows. This is currently still in an experimental state and we hope that widespread testing of 6.3.0 will mean the 6.4 release of WinGRASS will be fully functional and robust. Existing users will be happy to know that these new features do not disrupt the base GIS which remains as solid as ever and fully backwards compatible with earlier GRASS 6.0 and GRASS 6.2 releases.

Several infrastructure changes accompany this release with the project becoming a founding member of the Open Source Geospatial Foundation (OSGeo). This includes a new home for the website, the Wiki help system, source code repository, community add-on module repository, integrated bug tracking system, and formal membership for the project in a non-profit legal entity. We hope that these changes will guarantee that the GRASS community will be well supported and vibrant well into the future.

The Geographic Resources Analysis Support System (GRASS) is a Geographic Information System (GIS) used for spatial modeling, visualization of both raster and vector data, geospatial data management and analysis, processing of satellite and aerial imagery, and production of sophisticated presentation graphics and hardcopy maps. GRASS combines powerful raster, vector, and geospatial processing engines into a single integrated software package.

The GRASS GIS project is developed under the terms of the GNU General Public License (the GPL) by volunteers the world over. GRASS differs from many other GIS software packages used in the professional world in that it is developed and distributed by users for users, mostly on a volunteer basis, in the open, and is given away for free. Emphasis is placed on interoperability and unlimited access to data as well as on software flexibility and evolution rate. The source code is freely available allowing for immediate customization, examination of the underlying algorithms, addition of new features, and fast bug fixing.

GRASS is currently used around the world in academic and commercial settings as well as by many governmental agencies and environmental consulting companies.

Software download at and numerous mirror sites.

Full story at

GRASS GIS releases version 6.2.3

27 November 2007

The development team is happy to announce that a new bugfix version of GRASS GIS has been released today. This release fixes a number of bugs discovered in the 6.2.2 source code. It is primarily for stability purposes and adds minimal new features. Besides bug fixes it also includes a number of new message translations and updates for the help pages and projection database.

Highlights include further maturation of the GRASS 6 GUI, vector, and database code. Some improvements have been backported from the GRASS 6.3 development branch where new development continues at a strong pace of approximately one code commit every hour, including major work on an all new cross-platform wxPython GUI and a native MS Windows port (from 6.3.0 onwards).

The Geographic Resources Analysis Support System, commonly referred to as GRASS, is a Geographic Information System (GIS) combining powerful raster, vector, remote sensing and and geospatial processing engines into a single integrated software suite. GRASS includes tools for spatial modeling, visualization of raster and vector data, management and analysis of geospatial data, and the processing of satellite and aerial imagery. It also provides the capability to produce sophisticated 4D presentation graphics and hardcopy maps.

GRASS is currently used around the world in academic and commercial settings as well as by many governmental agencies and environmental consulting companies. It runs on a variety of popular hardware platforms and is Free open-source software released under the terms of the GNU General Public License.

GRASS is a proposed founding project of the new Open Source Geospatial Foundation. In support of the movement towards consolidation in the open source geospatial software world, GRASS is tightly integrated with the latest GDAL/OGR libraries. This enables access to an extensive range of raster and vector formats, including OGC-conformal Simple Features. GRASS also makes use of the highly regarded PROJ.4 software library with support for most known map projections and the easy definition of new and rare map projections via custom parameterization. Strong links are maintained with the QuantumGIS and R Statistics projects with integrated GRASS toolkits available for both.

Software download at and numerous
mirror sites.

Full story at