This is an old revision of the document!


Introduction to Open Source Geographic information systems: QGIS and GRASS

Quebec Centre for Biodiversity Science, http://qcbs.ca

Guillaume Larocque, Research Professional. (guillaume.larocque@mcgill.ca)

**THIS PAGE WAS DONE FOR AN OLDER VERSION OF THE WORKSHOP [[introgis|CLICK HERE FOR THE CURRENT VERSION]]**

Quantum GIS 1.8.0 and GRASS 6.4.2 (or 6.4.3) will be used for this workshop. Please make sure that you have the appropriate versions installed prior to the workshop.

For Windows users: Simply install QGIS and GRASS using the OSGeo4W installer available here.

For Mac users: Instructions to install QGIS and GRASS on Mac OS X are here. Please install GDAL and GRASS before installing QGIS.

For Mac OSX Leopard users: there are no complete packages available. The following packages have to be installed individually:

  * QGIS 1.8 "Leo"
  * GRASS 6.4.2
  * GDAL Framework 1.8.1-1 "Leo"
  * FreeType Framework 2.4.6
  * cairo Framework 1.10.2
  * UnixImageIO Framework 1.3.0
  * PROJ Framework 4.7.0
  * GEOS Framework 3.3.0
  * SQLite3 Framework 3.7.6.3

UIIO, PROJ, GEOS, SQLite3, and GDAL are available at: http://www.kyngchaos.com/software/frameworks#gdal_complete
All others available at: http://www.kyngchaos.com/software/archive

First, copy QGIS from the .dmg into your applications folder (don't try to run it yet). Now install the packages in the following order:
  * FreeType
  * PROJ
  * UnixImageIO
  * cairo
  * GEOS
  * SQLite3
  * GDAL
  * GRASS

For Linux users: Install QGIS and GRASS following these instructions.

General resources

Quantum GIS

GRASS

The Geographic Resources Analysis Support System (GRASS) is one of the oldest Geographic Information Systems still in existence today. It was originally developed in the early 1980s by the US Army as a software management tool for military applications. What was originally a command-line software running on Unix systems has now evolved into a multi-platform open source system with capabilities for raster and vector display and analysis, image processing, remote sensing, database management, and more. In recent years, GRASS became part of the Open Source Geospatial Foundation (OSGeo), a not-for-profit organization whose mission is to support and promote the collaborative development of open geospatial technologies and data. This foundation also supports multiple other projects for web mapping and desktop GIS, as well as various geospatial libraries for spatial data conversion, storage and interoperability. One of the projects supported by OSGeo is Quantum GIS (QGIS), a user-friendly, multi-platform, raster and vector desktop GIS which provides a point-and-click interface with functionalities for displaying, editing and analysing geospatial data. These functionalities can be greatly expanded through plugins that are either part of QGIS or that are supplied by third-parties. An interesting aspect of QGIS is the integration with GRASS, allowing the user to perform a number of GRASS operations from within QGIS. The combination of varied plugins and GRASS support turn QGIS into a full-fledge GIS with a strong analytical toolbox.

List of other Open Source or free GIS desktop programs of potential interest:

  • uDIG - An interesting and user friendly GIS written in Java for data access, editing, and viewing. It is under active development at the moment.
  • SAGA GIS - (System for Automated Geoscientific Analyses) A powerful GIS mostly aimed at data analysis and spatial algorithms.
  • ILWIS Open - A remote sensing and GIS software which integrates image, vector and thematic data in one package on the desktop.
  • gvSIG - A GIS currently in development having a user-friendly interface, being able to access the most common formats, and featuring a wide range of tools for query, layout creation, geoprocessing, networks, etc.
  • OpenJUMP GIS - The current version of OpenJUMP can read and write shapefiles other simple files, has limited support for the display of images and good support for showing data retrieved from web-services. It can be used as a GIS Data Viewer, or for editing geometry and attribute data.
  • Open QGIS - Click on Plugins­ > Manage Plugins, and activate fTools and GdalTools plugins.
  • Now you will see a dropdown menu in QGIS with options for vector analysis (fTools), raster analysis (GdalTools) and multiple other plugins (under plugins).
  • Go to Settings > Options ­­> CRS, and click on 'Enable on the fly reprojection by default'.
  • Start a New Project or close and reopen QGIS.

Layers in different reference systems will then be reprojected on the fly.

Click here to download a ZIP file containing all the files needed for the exercises below.

Exercise 1 - Digitizing and map display

Ojective: Create a map of the two McGill forests in the Montreal West Island: The Morgan Arboretum and the Molson Reserve. Draft a plan for the expansion of these reserves and create nice maps outlining your plan.

Step 1: Activate the OpenLayers Python plugin. If that isn't working for you, add this URL (name it OpenLayers) to the list of repositories : http://build.sourcepole.ch/qgis/plugins.xml , close the python plugins and try to enable the OpenLayers plugin 0.9.3 plugin.

Step 2: Open a text editor (for example, Notepad in Windows or TextEdit on Mac) and copy the following text to it.

ID,Name,Long,Lat
1,MorganArboretum,-73.9541,45.4303
2,MolsonReserve,-73.9763,45.3943

Then save as a text file with a .csv extension and add this as a vector layer to the QGIS map canvas with the >Layer>'Add Delimited Text Layer' function (if this isn't showing up for you, activate it in the plugins). Use Longitude as the X column and Latitude as the Y column. Specify >Geographic coordinate systems>WGS 84 as the CRS. Now, you should see two points on your screen indicating the center of each reserve.

Step 4: Add a Google Satellite layer using the Openlayers plugin (>Plugins>OpenLayers plugin>). Note that the CRS of the map canvas changes to Google Mercator (WGS 84 / Pseudo mercator) when the layer is added. Move this layer to the bottom of the maps canvas by clicking on it in the left menu and dragging it below the other layer.

Step 5: You need to add a New empty polygon layer (>Layer>New>Shapefile layer>Polygon). Specify the CRS as NAD 83/UTM zone 18N. A column for ID (integer) is already specified, add one for Name (Text data). Give the shapefile layer an appropriate name (ex. McGill_Reserves.shp). Enable Edit mode by left-clicking on the layer name in the left menu and selecting 'Toggle editing'. Digitize the McGill Molson Reserve and the Morgan Arboretum from the satellite image. To help you locate those areas, use the points which you imported from the CSV file. For the Molson Reserve, digitize the forested area located between the residential street (Blvd. Perrot) and highway 20 by clicking on the digitize icon: . For the Arboretum, digitize the contiguous forested area surrounding the point. The Arboretum is approximately 250 hectares. Now add polygons adjacent to those polygons showing what you propose as extensions to the existing reserves. For two adjacent polygons to properly share a common boundary, you need to digitize the second polygon so as to make it overlap slightly the first one. For this to work, it is important to select the 'Avoid int.' for that layer in the > Settings ­­­­> Snapping Options dialog.

Step 6: When done, exit Edit mode and save the layer.

Step 7: Start a new project, and set the CRS to NAD 83/UTM zone 18N.

Step 8: Add your new reserves shapefile to this project.

Step 9: From the Layers properties menu, change the colors of the Two reserves, add the names of the reserves as labels, and put appropriate labels for the extension of the reserves that you are proposing.

Step 10: Assign appropriate colors for each layer in the Style tab of the Layer Properties menu. Add the roads (Routes.shp), other forested areas (region_boise.shp) and water bodies (Region_hydrique.shp) to the map canvas. Then, using the Print composer (>File>New Print Composer), generate a map that shows the reserves you have digitized, complete with a title, a North Arrow, Labels, and a legend. Note: When you open the Print Composer, you are presented with an empty page. To add the map, you need to click on the Add New Map canvas and drag the area on the map where you want the map to appear.

CHALLENGE: Find an interesting WMS layer to add to your map.

Exercise 2 - Buffers and basic analyses

Objective: Find out the proportion of voters living within 500 m of a forested area larger than 10 hectares in the Jacques-Cartier provincial electoral district in the Montreal West Island.

Step 1: Start a new project and open the files named 'sections_vote_31h5.shp' and 'regions_boise.shp'. Convert these files to the NAD83/UTM 18N CRS by saving them as new files (right click on the layer name in the left menu>Save as).

Step 2: Start a new project again and set the project CRS to NAD83/UTM 18N. Add the two files you just saved.

Step 3: Open the Regions boisés attribute table and add a column containing the area in hectares of each forest. Toggle the layer editing and use the Field Calculator (small calculator icon at the bottom) to create the new column (use the Geometry>$area operator). Note that the UTM coordinate system is in meters and 1 hectare = 10,000 square meters. Important: Specify 'Decimal Number (real)', put the output width at 15 and increase the precision to 2.

Step 4: Perform a query (>Layer>Query) to isolate the forested areas larger than 10 hectares. Note that you can't perform a query while you are in Edit mode.

Step 5: Create a new shapefile with a 500m buffer surrounding the forested areas (>Vector>Geoprocessing tools>Buffer).

Step 6: Perform a query to select only the Jacques-Cartier electoral district (TRI_CEP column).

Step 7: Perform a 'Select by Location' query (in Vectors menu) to select all the voting zones that intersect the buffer.

Step 8: Use the 'Basic Statistics' function to find the sum of the 2008_12_08 column.

Step 9: Clear the selection and run the 'Basic Statistics' function again. Calculate the proportion by hand.

CHALLENGE : Find which of the electoral disctricts shown on the map contains the largest cumulative length of roads. Use the Dissolve and Sum line lengths (in Vectors menu) functions.

DAY 2

GRASS Terminology

  • DATABASE (Géodatabase): Main GRASS directory.
  • LOCATION (Secteur): Working directory. SCR unique.
  • MAPSET (Jeu de données): Sub-directory under the LOCATION.
  • REGION (Région): Définie par une étendue et une résolution spécifiques.
prefix Description Example
d.* display (graphical output) d.rast: views raster map, d.vect: views vector map
db.* database management db.select: selects value(s) from table
g.* general file operations g.rename: renames map
r.* raster data processing r.buffer: buffer around raster features, r.mapcalc: map algebra
i.* image processing i.smap: image classifier
v.* vector data processing v.overlay: vector map intersections
ps.* postscript map creation format ps.map: map creation in Postscript,
r3.* raster voxel data processing r3.mapcalc: volume map algebra

The full list of GRASS commands is here.

Key functions

g.list rast lists available raster maps
g.list -c vect lists available vector maps with attribute column names
g.proj -p Print projection information
r.info raster_name information about a raster
v.info vector_name information about a vector
g.region set current region
g.region rast=raster_name -p set region based on an existing raster

Exercice 3 - Interpolation and raster manipulation with GRASS, within QGIS

Objective: create two maps two compare the spatial distribution of ovenbird populations in Quebec from 1980-1994 and from 1995-2010, using the data from the Breeding bird survey (BBS).

IMPORTANT: Download this zip folder and move the files to your QGIS working folder, replacing the existing files you downloaded last week.

Step 1: Start a new project.

Step 2: Add the CSV file named 'BBS_Routes_QC.csv' to the map canvas using the function 'Add delimited text layer' (Choose 'Selected delimiters'→comma) and specify the CRS WGS84 (Geographic)

Step 3: Save this layer as a shapefile named 'BBS_Routes_QC.shp' with the CRS: NAD83 / Quebec Lambert.

Step 4: Add the file 'province.shp' to the canvas. Select the province of Quebec and save the polygon as a new shapefile with the CRS NAD83 / Quebec Lambert.

Step 5: Start a new project and specify the CRS NAD83 / Quebec Lambert. Add the files from steps 3 and 4 to the canvas.

Step 6: Add the (non-spatial) table 'BBS_Ovenbird_QC.csv' to the canvas using 'Layers>Add Vector Layer' and by listing all the file types (yes, this is counterintuitive…)

Step 7: From the menu 'Layer properties' of the table 'BBS_Routes_QC.shp', find the Join tab and join the table BBS_Routes and BBS_Ovenbird_QC using the shared column 'Route' (select Route as both the source and the target columns). You should now see the number of ovenbird observations per year in the attribute table of the 'BBS_Routes_QC' file.

Step 8: Activate the editing mode and use the 'Field Calculator' to add a new column (integer) to the table BBS_Routes containing the sum of the number of bird observations between 1980-1994 and 1995-2012. Each field has to be selected one by one from the Fields and Values list (sorry). The data type should be integer (length 10 or more). Column names should have less than 10 characters, should start with a letter and should contain no spaces or special characters

Step 9: Save the changes and exit editing mode. Note that only the new columns that were created are really part of the BBS_Routes_QC table. Now remove the join (from the Properties menu) to remove the yearly observations. If you see a series of NULL values in the newly computed columns, close the attribute table and open it again.

Step 10: Make sure that the GRASS extension is activated and that the GRASS icons bar is selected in >View>Toolbar, and that it is not hidden somewhere in your icons bar. Click on the New mapset icon and define a new GRASS Database, a new Location (name it workshop2), choose the CRS NAD83 / Quebec Lambert, specify the default GRASS Region by using the current QGIS extent and specify a name for the Mapset.

Step 11: You now need to import the files into GRASS. Click on 'Open GRASS tools', click on the 'Module tree'>'File management'>'Import into GRASS'>'Import a vector into GRASS'>'v.in.ogr.qgis'. Select the file BBS_Routes_QC.shp and specify a name for the output file. Add the newly created GRASS file to the map canvas by clicking on the 'Add Grass vector layer' icon.

Step 12: Click on 'Display current GRASS region' and on 'Edit current GRASS region'. You need to define the properties of the current GRASS region to obtain a resolution of 2km x 2km *exactly* so that you can include most of the points (you can exclude points in Northern Quebec), but without extending too much outside of the extent covered by the points. Don't hesitate to use a calculator and the image below:

Step 13: Convert the Quebec province file into GRASS format (same as Step 11).

Step 14: Find the v.to.rast.constant GRASS function and convert the Quebec polygon to raster (use 1 as Raster value). This will have the extents and resolution specified in Step 12.

Step 15: Use the r.mask GRASS function to create a mask with the raster created in Step 14. Add this new raster map to your canvas using 'Add GRASS raster layer'.

Step 16: Find the module v.surf.rst GRASS function, choose a name for the output file and interpolate a map for the data from 1980-1994 and another for the 1995-2005 data. You need to select the appropriate columns in the attribute field. Add these GRASS rasters to your map canvas. You will see that these maps will have the proper resolution (Step 12) and the interpolation will be limited by the mask (Step 15).

Step 17: In the properties of each layer, develop a new color palette appropriate for these files (Tab Style>Color palette>Color palette and then click on the tab >Colormap). Choose a linear color interpolation, define three relevant number values associated with each color (you can click on the raster with the 'Identify Features' tool to see the range of values) and assign the same palette to rasters from each period. Do you notice any difference in the spatial distribution of the ovenbirds over time?

CHALLENGE: Generate a scatter plot (in Excel or R) showing the relationship between the mean annual temperature (from the layer Quebec_mat_tenths.txt in degrees C x 10) and the ovenbird abundance at the BBS sites. You will need the extension 'Point sampling tool' which can be found in '>Vector>Analysis' once installed.

Download this zip file and make sure to extract the content in an easily accessible folder.

Exercise 4

Objective: create a NDVI (Normalized Difference Vegetation Index) map, associate it with a color palette and extract the mean NDVI for parcs and fields of the region.

Step 1: Within QGIS: Add the parcs_terrains_sports.shp layer to the canvas and save it as a new file with the CRS NAD83 / UTM 18N. Open the mapset 'landsat' from the 'Day 2' Database folder extracted from the zip file. Convert the file to GRASS format using the v.in.ogr.qgis function. If the operation was completed successfully, you can now close QGIS.

Step 2: Open GRASS. On the 'Welcome to GRASS' interface, Browse to select the GIS Data directory and select the folder named Day 2 that you have extracted from the downloaded zip file. The project should be set to '31h05' and the mapset to 'Landsat'.

Step 3: Visualize the Landsat bands by clicking on the 'Add raster map layer icon' on the 'GIS layer manager window' and selecting one of the landsat layers.

Step 4: Find the r.mapcalculator function under >Raster>Raster Map Calculator. Specify 'NDVI' as the name of the raster to create. You then need to find bands 3 and 4 in the 'Insert existing raster map' dropdown menu, and come up with this formula in the Expression box:

1.0*(band4-band3)/(band3+band4)

the 1.0 preceeding the formula assures that the resulting map will contain real numbers with decimals. Then click on Run. The NDVI raster should have been added automatically to the map canvas.

Step 5: Select >Raster>Manage colors>Color tables. Under the 'Required' tab, select the NDVI layer. In the Colors tab, select ndvi as Type of Color Table. The NDVI map should now change colors. On this map, dark green areas represent areas that are densely vegetated while blue areas are not vegetated.

Step 6: Add the layer containing the parcs and fields (parcs et terrains) by clicking on the icon 'Add vector layer'.

Step 7: Select Vector>Update attributes>Update attributes from raster. Choose the parcs layer as 'Name of vector polygon map' and the NDVI layer as 'Name of raster map. Use 'ndvi' as a column prefix.

Step 8: Click on the right mouse buton on the name of the layer and select 'Show attribute table'. By looking at the right part of the table, you will see that columns were added with the statistics of the NDVI cells for each field and parc.

Exercise 5

Objective: Create a 3D visualization showing a composite RGB map superimposed on an elevation surface.

Step 1: Select Raster>Manage colors>Create RGB. You will now create an RGB composite file by using the following band associations: Red:4 Green:3 Blue:2.

Step 2: In Files>NVIZ, create a 3D visualization with the file 31H05_dem as the elevation file and the RGB composite you just created (leave the other fieds blank). To improve the display resolution, click on Visualize>Raster Surface and specify 1 as the fine resolution.

*CHALLENGE* Create an interpolation elevation map from the Elevation_points.shp file and, using the v.surf.rst, function Compare this map with the 31H05_DEM file.

*CHALLENGE 2*: Visualize one of the vector files in R (basic instructions below).

Using R with GRASS

Before first use (in R):

install.packages('spgrass6', dependencies=TRUE)
install.packages('rgdal')

Mac users:
install.packages('spgrass6', type='source', dependencies=TRUE)
install.packages('rgdal', type='source', dependencies=TRUE)

spgrass 6 commands:

execGRASS               execute GRASS commands (when launched from R)
system                  execute GRASS commands (when launched from GRASS)
readRAST6               read GRASS raster files
writeRAST6              write GRASS raster files
readVECT6               read GRASS vector object files
writeVECT6              write GRASS vector object files
gmeta6                  read GRASS metadata from the current LOCATION
getLocationProj         return a PROJ.4 string of projection information
gmeta2grd               create a GridTopology object from the GRASS region
vInfo                   return vector geometry information
vColumns                return vector database columns information
vDataCount              return count of vector database rows
vect2neigh              return area neighbours with shared boundary length