NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Geospatial Data Introduction

Resources:

Understanding GRASS GIS commands

You can run GRASS several ways using its graphical user interface (GUI) and command line interface (CLI):

  • for display (d.*) commands use GUI, for other commands (r.*, g.*, v.*) type the command name into the command console in the Console tab (see the tabs at the bottom of the Layer Manager) and then press Enter to open command (module) GUI;
  • for the display (d.*) commands use GUI, for other commands type or paste (Paste or Shift Insert) command with all options into Console and Enter, you can use Up/Down arrows on your keyboard to retrieve previously run command, if you want to re-run it with different options;
  • run everything through GUI: start the command using a relevant button or select it from a drop down menu (use the Modules to find the command);
  • when running on Linux or Mac, paste the command with the options into the shell (you can copy/paste multiple commands at once)

Instructions for both GUI and command line approach are provided in the introductory assignments, but the instructions will be mostly command line later on. Feel free to use the approach that suits you the best. For a quick reference, here are the GUI equivalents for the d.* commands:

In GRASS GIS Layer Manager toolbar

In GRASS GIS Map Display toolbar The d.* commands throughout the assignments indicate which map layers should be included in the output map, d.out.file indicates which map should be saved for the report. Often, you will already have the required map layers in your Layer Manager, so you don't have to load them again using a d.* command. Note, that the CLI instructions can also be used for running the assignments on Linux and Mac as scripts.

GRASS project and mapset

If you don't already have GRASS GIS spatial database nc_spm_08_grass7 downloaded, you can find links to it on Course logistics page.

Before starting with GRASS GIS it is important to understand that GRASS GIS uses the concept of projects, mapsets and computational region to support efficient analysis and modeling without the need to subset and resample data. Data in a GRASS project have common coordinate system. Each project can have several mapsets, one is called PERMANENT to store baseline data, others are set by users. You can modify the data only in your current mapset, but you can view or copy data from other mapsets. See also structure of GRASS database.

DO NOT MOVE the files under project directory using non-GRASS tools.

Suggested workflow is to create a new mapset for each assignment and do all computations there.

Computational region is set by the command g.region and specifies the spatial extent and resolution at which raster operations are performed.

Start GRASS GIS

Start GRASS GIS as any other program on your computer. You should see a project nc_spm_08_grass7 in your database. If you have not downloaded it yet, click the icon Download sample project to current database and select Complete North Carolina dataset, click Download, close when finished.

To start working on your Assignment 1 right click on nc_spm_08_grass7 and select Create mapset, type in name for the new mapset, for example, assignment1, press OK. You will get a message that the Current mapset is assignment1 and any new map layer that you will create will be saved here. The provided data are in the mapset PERMANENT - double click to see the list of available raster and vector map layers. You can preview a map layer by double-clicking on it but we will manage the Map Display in the Layers tab and/or using display commands.

To save your outputs and to store external data create a directory where you can write, e.g. on MS Windows, C:\Users\myname\Documents\GIS582_HW\report1. Then set this directory in GRASS GIS as your working directory:
Settings > GRASS working environment > Change working directory > select this directory
or type cd (stands for change directory) into the GUI Console and hit Enter: 

cd

Find the working coordinate reference system, spatial extent and resolution:
In GUI: Settings > Map projections > Display map projection
In GUI: Settings > Region > Display Region

Or paste the following command into the Console tab: 

g.proj -p
g.region -pl

Display data in 2D

First set working region to a given elevation raster using GUI:
Settings > Region > Set Region > Set region to match raster map(s) > select elevation > Run > Close or alternatively, paste in the command console in the Console tab: 
g.region raster=elevation -p

Display elevation map:

In Layer Manager > Add raster layer button > select raster "elevation" > OK.
In Map Display > Zoom to computational region extent.
Make sure the Render button (lower right corner) for automatic rendering of maps is checked on. You can also use the first button Display map in Map Display to render maps.

Alternatively, you can add raster layer by pasting this command into the Console tab: 

d.rast elevation

Display elevation as a colored shaded map:

In Layer Manager > Add various raster maps layers > Add shaded relief map layer > select "elevation_shade" as shaded relief map and "elevation" as color map.
Under Optional tab set Percent to brighten to 30 > OK.

Command line equivalent: paste into the Console and press Enter: 

d.shade shade=elevation_shade color=elevation brighten=30

Now check the range of elevation values of the entire map. Either by pasting the following command into the Console or right clicking on the elevation raster layer and selecting Metadata in the context menu. 

r.info -r elevation

Display vector line and point maps:

In Layer Manager: Add vector map layer button > select "streets_wake" > OK. Adjust the colors, line widths, symbols by clicking on relevant tabs. Add additional maps ("roadsmajor", "schools_wake").

These are the command line equivalents: 

d.vect streets_wake
d.vect roadsmajor color=red width=2
d.vect schools_wake icon=basic/circle size=10 fill_color=blue

Add legend and scale bar:

To add the legend use Add map elements button in Map Display and select Add legend (you can use mouse to move it around, right click to adjust the size or remove it).a
To add scale bar use: Add map elements > Add scale bar (double click on it to change it, for example you can change its length and units under Optional tab).

These are the command line equivalents: 

d.legend elevation
d.legend.vect
d.barscale

Save the displayed map:
In Map Display click Save display to graphic file. Or run the equivalent command from the console, the graphic file will be saved into your working directory. 

d.out.file mymap

Display data in 3D

This is just a quick preview: We will explore 3D visualization in future assignments.

In Layer Manager check off the elevation_shade layer, but keep the elevation raster on. In Map Display switch 2D view to 3D view. After the 3D mode is loaded, increase viewing height with slider, rotate puck to the south. Click Data tab and set Fine mode resolution to 1, to get the full resolution surface. Save the image by clicking on Save display to graphic file (tiff support is required). You may need to install IrfanView to view TIFF images saved from the 3D view. Alternatively, you can just take screenshots of the 3D visualization.

Changing data view, query and region

To get familiar with query and zoom tools, set the view back to 2D. Querying data means finding out the values in raster maps or attributes in vector maps at a location selected by mouse. In Layer Manager select the map layer(s) to query by clicking on it. In Map Display click button Query raster/vector map(s) and click on a map at locations where you want to know the values/attributes.

In Map Display explore the options to zoom: zoom-in, zoom-out, zoom to selected map, zoom to computational region. You can also zoom to specific map layer by right clicking on the layer and selecting Zoom to selected map(s).

Next, display statewide DEM, county boundaries and climate stations:

Remove all map layers from Layer Manager, remove legend and scale bar by right clicking on it and select Remove.
In Layer Manager menu: Settings > Region > Set Region.
In g.region dialog: Set current region from named region > select "nc_500m" > Run > Close.
In Map Display: Zoom to computational region.
In Layer Manager: Add raster layer > select "elev_state_500m".
Add vector map layer "precip_30ynormals", set the symbol to basic/marker.
Add vector map layer boundary_county, to display only county boundary, under tab Selection switch off area and switch on boundary. Rearrange the layers by mouse if needed and add legend and scale bar as in the previous example.

The command line equivalent is belowr:, 

d.erase
g.region nc_500m -p
d.rast elev_state_500m
d.vect precip_30ynormals icon=basic/marker
d.vect boundary_county type=boundary
d.out.file nc_precip_stations

Note that any new data (map layers) that you will import in the next section or create in the future assignments are automatically stored in your GRASS database/project/current_mapset. If you want to go back to displayed maps that you have worked on in Map Display, you can save the Map Display and Layer manager, including the symbology and layout in a GRASS workspace file *.gxw under File > Workspace > Save as.

Import and export for raster and vector data

In this section we will cover importing and exporting of raster and vector data which are in the same projection as the GRASS project we work in. Cases where the projection differs are covered in a separate guide.

Import and export Shapefiles

Download a shapefile with geodetic points: geod_pts_spm.zip. Save and extract (unzip) it to your working directory.

To import file in GUI use:

  • File > Import vector data > Simplified vector import with reprojection.
  • Use Browse button to find the "gdc.shp" file you extracted earlier. to fill the Source input field.
  • Then click on Import button.
In case the imported point data was not displayed automatically, use File > Add vector to display the data and zoom to the vector layer.
How are the imported geodetic data points spatially related to the "streets_wake" vector data? Hint: add "streets_wake" to layer manager and zoom to this layer.

Export the roadsmajor map as a shape file:

  • File > Export vector map > Common export formats
  • Select "roadsmajor" as Name of input vector map
  • Use Browse to find a directory to export the data to, e.g. your "working directory". This will fill out the Name of OGR output datasource field with the path to this directory (e.g. "C:\mydirectory\roadsmajor.shp") Make sure to specify the full path to your output files otherwise you might have hard time finding them especially on MS Windows. (Also note that you need to have write permissions for that directory.)
  • Select ESRI_Shapefile in Data format to write
  • click Run
You can then display the exported "roadsmajor" data in ArcGIS. Note that the names of the GRASS GIS tool used in the background during the import and export were: v.import and v.out.ogr.

Import a raster file provided in GeoTIFF format

Download a landuse raster in GeoTIFF format: lc96ras_cut.tif and save it into your working directory.

To import this file through GUI use: File > Import raster data > Simplified raster import with reprojection. Use Browse button to find the path to the downloaded "lc96ras_cut.tif" file.
Make sure the file is checked in the List of GDAL layers. Click Import.

Alternatively, if the data that you want to import is stored in your working directory you can use the following command: 

r.import input=lc96ras_cut.tif output=landuse96_subset
Note, that if the data is not in the working directory you need to provide a full path to your input data, so you would have to replace "lc96ras_cut.tif" by something like "C:\path_to\my_file\lc96ras_cut.tif". In such case the import is easier to do through GUI where you can use Browse to get the correct path to your data.

Exporting raster files as ASCII grid and GeoTIFF

Now, export raster files as ASCII grid and GeoTIFF.

Use File > Export raster map > Common export formats.
Similarly as for the vector data export, the exported file will be savedi in your working directory or you need to specify the full path to the directory where you want the exported file to be saved. Use GTiff and AAIGrid as Raster data format to write and use file extension appropriate for each format, i.e., .tif and .ascii respectively.
Alternatively, you can use the command line equivalent to export the data to your working directory: 

r.out.gdal input=basin_50K output=basin_50K.tif
r.out.gdal input=elev_ned_30m output=elev_ned_30m.ascii format=AAIGrid

Creating new GRASS project

We will create new projects based on EPSG code and based on georeferenced file.

First we create a new project with Coordinate Reference System (CRS) given by EPSG code 3404. Which CRS has EPSG 3404 and what are its characteristics?

Launch Project Wizard from Layer Manager menu > Settings > GRASS working environment > Create new project
In Name field type "nc_spf" > Next
Select method Select CRS from a list by EPSG or description > Next
Type EPSG code 3404 > Next > OK > Finish
In this new project review the CRS information: 

g.proj -p

Another way to set up a GRASS project uses CRS information from a georeferenced file. Again launch the Project Wizard:

In Name field type "nc_spm" > Next
Select method Read CRS from a georeferenced data file > Next
Browse to downloaded file lc96ras_cut.tif > Next > Finish
Import data? > Yes > OK
Review again the projection information: 

g.proj -p

Optional: add WMS layer

When we created new projects in the previous section (nc_sft and nc_spm) GRASS automatically switched to the mapset in a newly created project (it will be bold in the data catalog). For the next task we go back to our "assignment1" mapset in the "nc_spm_08_grass7" project:

In the Data catalog, Right click on the "assignment1" in the "nc_spm_08_grass7" project and select Switch mapset from the menu. You should now be ready to set the region to a small rural watershed and display an existing ortho and then view and download orthophoto through web mapping service. Requires good internet connection.

First change region to the small rural area specified by the saved region named "rural_1m", display the provided orthophoto and save its image for your report. 

g.region region=rural_1m -p
d.rast ortho_2001_t792_1m
d.out.file ortho2001

In Layer Manager toolbar Add web service layer.
Paste the link to the service into server field and press connect:

https://imagery.nationalmap.gov/arcgis/services/USGSNAIPPlus/ImageServer/WMSServer?
Wait until GUI changes. Be patient, it can take up to one minute.
In Available web services, select WMS 1.3.0.
In List of layers, there should be "USGSNAIPPlus", click on plus icon to open list of provided products.
Select a desired product (e.g. Natural colors or NDVI), Press Add layer button and wait.
The dialog can be closed afterwards.

To save the layer:
Right click on the web service layer in Layer Manager - choose Save web service layer. In dialog, set the name of the layer ortho_new and press Save layer.
Save the image of the new orthophoto for your report: 

d.out.file orthonew