Introduction to Open Source Geographic information systems: QGIS and GRASS

Quebec Centre for Biodiversity Science.

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

November 2011

Quantum GIS 1.8.0 and GRASS 6.4.2 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 standalone installer available here. 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.
  • Click on Plugins > Fetch Python Plugins­­­ > Repositories and click on 'Add 3rd party repositories'.

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 maps 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.

Step 2: Create a text file with a .csv extension in a text editor and copy the following text to it.

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

Then, add this as a vector layer to the canvas with the 'Add Delimited Text Layer' function. Use Longitude as the X column and Latitude as the Y column.

Step 3: Add a Google Satellite layer using the Openlayers plugin. Notice that the display switches to the Google Mercator projection. Move this layer to the bottom of the maps canvas.

Step 4: Digitize the McGill Molson Reserve and the Morgan Arboretum from the satellite image. You need to add a New polygon layer, enable edit mode, add a column for ID and Name. 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. 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. It is important to select the 'Avoid int.' for that layer in the Settings ­­­­> Snapping Options dialog.

Step 5: Exit Edit mode and save the layer in the UTM zone 18N reference system.

Step 6: Start a new project, and set the reference system of the project (CRS) to UTM zone 18N.

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

Step 8: 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 9: Generate a map that shows the reserves you have digitized, roads (Routes.shp), other forested areas (region_boise.shp) and water bodies (Region_hydrique.shp). Assign appropriate colors for each layer in the Style tab of the Layer Properties menu. Make a nice looking map using the print composer, complete with a title, a North Arrow, Labels, and a legend.

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

Exercise 2 - Buffers and basic analyses

Objective: Find out how what 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.

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: Add a column to the Regions boisés attribute table containing the area in hectares of each forest (note: 1 hectare = 10,000 square meters).

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 (choose only the selected features).

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

Step 7: Perform a Select by Location (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 Join Attributes by Location (in Vectors menu) function.

Terminologie GRASS

  • 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

—-



DAY 2

DOWNLOAD THE ZIP FILE needed for Exercises 3 and 4.

Exercise 3 - Interpolation and basic raster operations

Objective: create two maps comparing the spatial distribution of ovenbird populations in Quebec from Breeding Bird Survey route data from 1980-1995 and 1995-2010. Step 1: Start a New Project.

Step 2: Load the CSV file named 'BBS_Routes_QC.csv' into the canvas using 'Add Delimited Text Layer' (comma delimited).

Step 3: Load the 'province.shp' file into the canvas.

Step 4: Select the province of Quebec polygon as a separate shapefile in the Quebec Lambert (NAD83) CRS. Also save the BBS_Routes_QC file as a shapefile in the same CRS.

Step 5: Start a New project, specify the project CRS as Quebec Lambert (NAD83). Load the two shapefiles from step 4.

Step 6: Load the table named BBS_Ovenbird_QC.csv into the canvas using 'Add vector layer' and listing all file types.

Step 7: From the Layer Properties dialog, find the Join tab and join the BBS_Routes table with the BBS_Ovenbird_QC tables using the shared Route column.

Step 8: Use the Field Calculator to add a new column to the BBS_Routes table that contains the summed number of birds seen between 1980-1995 and 1995-2010.

Step 9: Once the new fields are created, save your edits, exit edit mode.

Step 10: Use the Grid (interpolation) function under Raster ­> Analysis to create an interpolated raster map with the Ovenbird distribution for each time period. Use the default interpolation options and define the extent of the map so that it encompasses all the BBS routes. Change the values for Extent and Size to have a raster that has 2 km pixels. Use the same properties for each raster.

Step 11: Clip each raster layer using the Clipper function and using the Province of Quebec shapefile as a mask layer.

Step 11: In the layer properties choose Colormap as a Style (Style tab) and generate your own three-color colormap from the Colormap tab, by using the raster values as an indication of the values that should be assigned to each color. Make sure to use the Linear Color interpolation option. Specify the same colormap for the two maps.

CHALLENGE: Generate a scatterplot (in Excel or R) showing the relationship between Mean Annual Temperature (from the Quebec_mat_tenths.txt layer in degrees C x 10) and ovenbird abundance, at the BBS routes locations from 1980-2010. Hint: you might need the 'point sampling tool' plugin.

GRASS Terminology

  • DATABASE: Main GRASS working folder.
  • LOCATION: Working directory. Unique map projection.
  • MAPSET: Subdirectories under location.
  • REGION: Defined by a specific extent and resolution.
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

Find a Full list of GRASS commands here

A few 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

Exercise 4 - Basic Raster Analysis with GRASS

Objective: Create a 3D visualization of an elevation map created with elevation values from ground control points, excluding water bodies. Learn how to import layers into R, perform basic manipulations in GRASS within QGIS and GRASS standalone.

Step 1: If you haven't done so already, save the 'region_hydrique' in the UTM 18N (NAD83) CRS. Start a New project and set the CRS of the current canvas to UTM 18N (NAD83). Add the region_hydrique file (in UTM) and the 'Elevation_points.shp' file in QGIS (this file is already in UTM).

Step 2: Make sure that the GRASS plugin is activated and that the GRASS toolbar is selected under View­­ > Toolbar. Click on the 'New Mapset' icon and define a new GRASS database folder, create new location (give it a name such as 'Workshop2', choose the UTM 18N (NAD83) CRS, set the region to the current QGIS extent, and choose a name for the mapset (e.g. 31h05).

Step 3: You now have to import the files into GRASS. Click on the 'Open GRASS tools icon', click on 'module tree', and find File management > import into GRASS > import vector into GRASS > v.in.ogr.qgis - Import loaded vector. Select the vector layer region hydrique (in UTM) and find name for the output vector map. Repeat these steps for the Elevation_points vector file. Load the newly created GRASS files into the canvas by clicking on the 'Add GRASS vector layer' icon.

Step 4: Click on the 'Display current GRASS region' icon and then on the 'Edit Current GRASS region' icon, and select a GRASS region that just encompasses all the elevation points. Set the Cell width and height to 30 m.

Step 5: Open GRASS tools and find the v.to.rast.constant module to convert the water bodies file to raster. Keep the raster value at 1.

Step 6: Close QGIS. We could decide to continue working within QGIS, but instead, we will move to GRASS outside of QGIS.

Step 7: Open GRASS and select the Database, Location and Mapset that you have created within QGIS.

Step 8: Load the Raster map you have created in Step 5, and the elevation points you have imported in Step 3 with the 'Add raster map layer' and 'Add vector layer' icons.

Step 9: Use the Raster ­­­> Mask module (r.mask) to create a mask layer for the next step. Select the 'Create inverse mask' and put 1 in 'Category Value to use for Mask'. Click the 'Run' button to run the module.

Step 10: Use the v.surf.rst module (Raster > interpolate surfaces) to create an interpolated elevation map using the Elevation_points layer. Specify a name for the output elevation map in the output tab and specify the ELEVATION column as the 'Name of the attribute column to be used for approximation'.

Step 11: Generate a slope map from this new elevation layer using the r.slope.aspect module (Raster > Terrain Analysis). Specify an output name for the slope layer.

Step 12: View the resulting map in three dimensions in the NVIZ module (File > NVIZ).

Using R with GRASS

On first use (within R):

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

For Mac users, run instead:
install.packages('spgrass6', 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