Geospatial Modeling and Analysis

Data display and visualization



In startup pannel set GRASS GIS Database Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata. For GRASS Location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for GRASS Mapset create a new mapset (called e.g. HW_data_visualization) and click Start GRASS session.

Change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter:


Changing the default font

Change the default font used for map rendering, use: Settings > Preferences > Map Display. Pick a new default font and Save the settings.

Visualization in 3D perspective

Set the region to the raster "elevation".
g.region raster=elevation -p 

Interactively view "elevation" raster and vector data ("streams", "usgsgages") in 3D view, adjust viewing position, surface properties and lighting to highlight features. Note, that legend is not currently supported in 3D view.

Remove or switch off any map layers in the Layer Manager. Then add raster "elevation", Add vector "streams", Add vector "usgsgages".
Zoom to computational region and Switch to 3D view.

Follow the video Visualization I and II for similar tasks (link is above) and save 2-3 images for your report (save as tiff or take screenshot, if you don't have tiff support).

Visualizing multiple surfaces

Visualize multiple surfaces (bare earth and surface with vegetation and structures),
analyze their relationship using crossections generated by interactive cutting planes

It is recommended to quit GRASS before starting the task below because we will be working with a smaller, high resolution region.

Start GRASS7 with your previous MAPSET. First set region to "rural_1m", then interpolate surface with vegetation from multiple return lidar points using the module (we will explain interpolation later on).

g.region rural_1m -p input=elev_lidrural_mrpts elevation=elev_lidruralmr_1m npmin=120 segmax=25 tension=40 smooth=1.0

In Layer Manager Add rasters "elev_lidruralmr_1m" and "elev_lid792_1m" (bare earth). In Map Display Zoom to computational region.

Make sure that you have only the two elevation maps switched on in the Layer Manager. Switch to 3D view and follow the video Visualization III (link is above). Save 2-3 interesting images for your report, include at least one crossection.

Basic 2D display operations

GUI is recommended for the tasks below, see the GUI equivalents for selected d* commands above, the command line is to indicate the workflow and output.

Display subsets of data

Visualy explore relation between developed areas and topography.

Display land use categories 1, 2 (developed land) over shaded topography.
Settings > Region > Set region > Set region to match this raster > select "landuse96_28m".
Add raster > select "elevation_shade".
Add raster > select "landuse96_28m" > Selection > List of cats to display > 1,2.
Zoom to computational region.
Right click on "landuse96_28m" map layer and Change opacity to show topography blended with landuse. Save display to graphic file. You may want to remove all previously used layers from the Layer Manager before you start visualizing the new ones.

g.region raster=landuse96_28m -p
d.rast elevation_shade
d.rast landuse96_28m values=1,2
d.out.file landuse_elev

Change colors for raster maps

There are many ways how to do it see r.colors manual).

First create your own copy of the elevation map (see the command below).
Display it by Add raster > "myelev".

g.copy raster=elevation,myelev
d.rast myelev

Set the color table to a predefined one, for example "bgyr" (blue, green, yellow, red).

GUI options: Right Click on the raster layer "myelev" in Layer Manager > Set color table > Define > Name of color table > "bgyr" > Run.
If you don't see the new colors - click the second button in Map Display to redraw.

r.colors myelev color=bgyr
d.rast myelev

Change the color table interactively using rules:

In the Layer Manager manu click Raster > Manage colors > Manage Color rules interactively.
Select raster map "myelev".
Click on colored squares to change the color, click Preview.
If it looks good - click OK.

To create color table for future multiple uses type the rules into a plain text mycolor.txt file using a text editor (TextEdit, Notepad), for example:

50 blue
70 aqua
90 green
110 yellow
130 orange
160 brown

Then assign it as follows:
Right click on "myelev" map layer > set color table > Define > Path to rules file > Run.
If needed, redraw using the second button on Map Display. You can also create and assign the above color table by typing or pasting it into the enter values interactively window or run it from comand line (make sure you have the correct path, in the following example mycolor.txt should be replaced by e.g. /Users/john/mycolor.txt on Mac OS X):

r.colors myelev rules=mycolor.txt

Note that you can combine color RGB and name definitions as well as % and cat/val to create complex custom color tables and store them in a text file for future use, see the examples in r.colors man page.

Compare the use of equal interval and histogram equalized color table for slope

To add the legend use Add map elements button on Map Display.
g.copy raster=slope,myslope
r.colors myslope color=bgyr
d.rast myslope
d.legend myslope
d.out.file myslopecolor

To explain the difference between the two maps, you can generate a histogram. On map Display click on Analyze map > Create histogram to open the histogram tool and save results to graphics file. Or use this command which will add a histogram as a layer to Layer Manager.

d.histogram myslope

Now set the histogram equalized color table, and save the new slope map and histogram images. Describe the effect of the histogram equalized color table.

r.colors -e myslope color=bgyr
d.out.file myslopecolorequalized
d.histogram myslope

Modify legend, scale and grid

To re-size the legend for myslope, right click on the legend, select Resize legend and resize with mouse. Alternatively you can resize by selecting Pointer mode in Map Display, double clicking on legend to launch the legend dialog and set options > Optional > Placement > 50,90,4,7.
Numbers are bottom,top,left,right as percentage of screen coordinates.
Add units to the legend: Add map elements > Add text layer > type deg > OK.
Add barscale on Map Display: Add map elements > Show/hide north arrow (double click on it to change it).

Note: you can use horizontal legends by using Placement at=6,10,2,30 or just stretch it horizontally with mouse.

d.legend myslope at=50,90,4,7

Switch off "myslope", display "myelev" raster and change the legend to "myelev" raster. Add grid for state plane coordinates at 5000m with ticks at 1000m. Also add a lat/long grid at 2 arc minute interval.

In Layer Manager > Add various overlays > Add grid layer > Size 5000 > Run.
To draw only border with ticks: change Size to 1000 and click Disable > Disable grid drawing > Run.
Switch off Disable grid drawing.
To draw lat/long grid: set Required > Size to 0:02 > click Draw > Draw geographic grid > Run.
Command line alternative:

d.rast myelev
d.legend myelev
d.grid size=5000 color=brown
d.grid -n size=1000
d.grid -g size=0:02 color=black
d.out.file myelevmap

Display color composite

In Layer Manager: Add various raster map layers > Add RGBmap layer.
d.rgb red=lsat7_2002_30 green=lsat7_2002_20 blue=lsat7_2002_10
d.vect roadsmajor color=yellow
d.out.file mylandsat

Export NC precipitation stations data to KML and post in Google Earth

To see what you are exporting, add the precip_30ynormals raster or use the following command:
d.vect precip_30ynormals
You may need to use right click in Layers to zoom to the extent of the whole vector map. Then use v.out.ogr to export the vector map in the KML format. If your are running the command in GUI do not forget to click point for data type. To view the data in Google Earth just drag and drop the KML file in the Google Earth window.
Make sure you include extension ".kml" in the output file name.
v.out.ogr precip_30ynormals output=precipitation.kml format=KML type=point
You need to set the file path to the desired directory or work using the current working directory. On MS Windows, you may do something like this:
v.out.ogr precip_30ynormals output=C:\temp\precipitation.kml format=KML type=point

Optional part

Add WMS layer

View and download orthophoto through web mapping service. Requires good internet connection. If it doesn't work, report the problem including any error message in the report.
First change region to the small rural area specified by the saved region named "rural_1m" and display the provided orthophoto.
g.region region=rural_1m -p
d.rast ortho_2001_t792_1m

In Layer Manager toolbar Add web service layer.
Paste the link to the service into server field and press connect:
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 "USGSNAIPImagery", click on that.
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.

Map Swipe

Use mapswipe to compare the new orthophoto with the ortho provided in the data set:
File > Map Swipe.
Add raster > select "ortho_2001_t792_1m" and "ortho_new".

Move slider to compare images - what changes do you observe? Save a screenshot of the mapswipe window showing the comparison of orthos. See more options here.