Difference between revisions of "Analytical data visualizations at ICC 2017"

From GRASS-Wiki
Jump to navigation Jump to search
Line 256: Line 256:
[[File:Elevation_slope_3D_view.png|800px|thumb|center| 3D visualization of elev_lid792_1m DEM with slope draped over]]
[[File:Elevation_slope_3D_view.png|800px|thumb|center| 3D visualization of elev_lid792_1m DEM with slope draped over]]


== Scripting in GRASS GIS ==
== Introduction to scripting in GRASS GIS with Python ==
 
The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python editor integrated in GRASS GIS accessible from the toolbar or the Python tab in the Layer Manager. Another option is to write the Python code in your favorite plain text editor like Notepad++ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu File -> Launch script.
 
 
<center><gallery perrow=3 widths=380 heights=190>
Image:Simple python editor v buffer.png|''Simple Python Editor'' integrated in GRASS GIS (since version 7.2) with ''Python'' tab in the background which contains an interactive Python shell.
Image:GRASS GUI Python shell.png|''Python'' tab with an interactive Python shell
</gallery></center>
 
 
The GRASS GIS 7 Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include:
* '''{{pyapi|script|script.core|run_command}}''': most often used with modules which output raster/vector data where text output is not expected
* '''{{pyapi|script|script.core|read_command}}''': used when we are interested in text output
* '''{{pyapi|script|script.core|parse_command}}''': used with modules producing text output as key=value pair
* '''{{pyapi|script|script.core|write_command}}''': for modules expecting text input from either standard input or file
 
Besides, this library provides several wrapper functions for often called modules.
 
==== Calling GRASS GIS modules ====
We will use GRASS GUI Simple Python Editor to run the commands. You can open it from Python tab.
For longer scripts, you can create a text file, save it into your current working directory and run it with <tt>python myscript.py</tt> from the GUI command console or terminal.
 
When you open Simple Python Editor, you find a short code snippet.
It starts with importing GRASS GIS Python Scripting Library:
<source lang="python">
import grass.script as gscript
</source>
 
In the main function we call {{cmd|g.region}} to see the current computational region settings:
 
<source lang="python">
gscript.run_command('g.region', flags='p')
</source>
 
Note that the syntax is similar to bash syntax (<code>g.region -p</code>), just the flag is specified in a parameter. Now we can run the script by pressing the Run button in the toolbar. In Layer Manager we get the output of g.region.
 
Before running any GRASS raster modules, you need to set the computational region. In this example, we set the computational extent and resolution to the raster layer <tt>elevation</tt>.
Replace the previous g.region command with the following line:
 
<source lang="python">
gscript.run_command('g.region', raster='elevation')
</source>
 
The run_command() function is the most commonly used one. We will use it to compute viewshed using {{cmd|r.viewshed}}.
Add the following line after g.region command:
 
<source lang="python">
gscript.run_command('r.viewshed', input='elevation', output='viewshed_python', coordinates=TODO, overwrite=True)
</source>
 
Parameter <code>overwrite</code> is needed if we rerun the script and rewrite the raster. Now let's look at how big the viewshed is by using
{{cmd|r.univar}}. Here we use <code>parse_command</code> to obtain the statistics as a Python dictionary
 
<source lang="python">
univar = gscript.parse_command('r.univar', map='viewshed_python', flags='g')
print univar['n']
</source>
 
The printed result is the number of cells of the viewshed.
 
List of convenient wrapper functions with examples:
* Raster metadata using {{pyapi|script|script.raster|raster_info}}: <code>gscript.raster_info('viewshed_python')</code>
* Vector metadata using {{pyapi|script|script.vector|vector_info}}: <code>gscript.vector_info('roads')</code>
* List raster data in current location using {{pyapi|script|script.core|list_grouped}}: <code>gscript.list_grouped(type=['raster'])</code>
* Get current computational region using {{pyapi|script|script.core|region}}: <code>gscript.region()</code>


== Introduction to Blender ==
== Introduction to Blender ==

Revision as of 02:53, 20 June 2017

Icc2017 logo.png
Grassgis logo colorlogo text whitebg.png

This is material for two sessions at the ICC 2017 workshop called Analytical data visualizations with GRASS GIS and Blender and Mapping open data with open source geospatial tools held in Washington, DC, July 1-2, 2017. These two sessions introduce GRASS GIS, example of its processing capabilities, and visualization techniques relevant to spatio-temporal data and high performance computing (HPC). Participants interactively visualize open data and design maps with several open source geospatial tools including Tangible Landscape system.

Software

We will use GRASS GIS 7.2.

MS Windows

Then download standalone GRASS GIS 7.2. binaries (64-bit version, or 32-bit version) from grass.osgeo.org. During the installation you can also download North Carolina sample dataset so please select this option. You can also download the dataset later (see the following section).

Mac OS

Install GRASS GIS 7.2 using Homebrew osgeo4mac:

brew tap osgeo/osgeo4mac
brew install grass7

OSGeo-Live

All needed software is included except for Blender.

Ubuntu

Install GRASS GIS from packages:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install grass

Linux

Install GRASS GIS from packages:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install grass

For other Linux distributions other then Ubuntu, please try to find GRASS GIS in their package managers.

Data

Prepared data

  • digital surface model derived from North Carolina Flood Plain Mapping lidar
  • orthophoto

Download OpenStreetMap data

We will use overpass turbo (web-based data filtering tool) to create and run Overpass API queries to obtain OSM data.

You can use Wizard to create simple queries on the currently zoomed extent. Zoom to a small area and paste for example this into the Wizard and run the query:

highway=* and type:way

The query was built in the editor and ran and you can see the results in a map now.

Now, paste this query in the editor window:

[out:json][timeout:25];
// gather results
(
  // query part for: “highway=*”
  way["highway"](35.76599,-78.66249,35.77230,-78.65261);
);
// print results
out body;
>;
out skel qt;

It finds roads in our study area. When you run the query, the roads appear on the map and we can export them as GeoJSON (Export - Data - as GeoJSON).

Now let's also download all building footprints and again export them as GeoJSON:

[out:json][timeout:25];
// gather results
(
  // query part for: “building=*”
  way["building"](35.76599,-78.66249,35.77260,-78.65261);

);
// print results
out body;
>;
out skel qt;

In the next steps we will import this data into GRASS GIS.

Introduction to GRASS GIS

Here we provide an overview of GRASS GIS. For this exercise it's not necessary to have a full understanding of how to use GRASS GIS. However, you will need to know how to place your data in the correct GRASS GIS database directory, as well as some basic GRASS functionality.

Structure of the GRASS GIS Spatial Database

GRASS uses unique database terminology and structure (GRASS GIS Spatial Database) that are important to understand for working in GRASS GIS efficiently. You will create a new location and import the required data into that location. In the following we review important terminology and give step by step directions on how to download and place you data in the correct place.

  • A GRASS GIS Spatial Database (GRASS database) consists of directory with specific Locations (projects) where data (data layers/maps) are stored.
  • Location is a directory with data related to one geographic location or a project. All data within one Location has the same coordinate reference system.
  • Mapset is a collection of maps within Location, containing data related to a specific task, user or a smaller project.

Creating a GRASS database for the tutorial

Unless you already have a directory called grassdata in your Documents directory (on MS Windows) or in your home directory (on Linux), create one. You can use the Browse button and the dialog in the GRASS GIS start up screen to do that.

Displaying and exploring data

Now that we have the data in the correct GRASS database, we can launch the Graphical User Interface (GUI) in Mapset user1.

The GUI interface allows you to display raster, vector data as well as navigate through zooming in and out. More advanced exploration and visualization is also possible using, e.g., queries and adding legend. The screenshots below depicts how you can add different map layers (left) and display the metadata of your data layers.

GRASS GIS modules

One of the advantages of GRASS is the diversity and number of modules that let you analyze all manner of spatial and temporal. GRASS GIS has over 500 different modules in the core distribution and over 230 addon modules that can be used to prepare and analyze data layers.

GRASS functionality is available through modules (tools, functions). Modules respect the following naming conventions:

Prefix Function Example
r.* raster processing r.mapcalc: map algebra
v.* vector processing v.clean: topological cleaning
i.* imagery processing i.segment: object recognition
db.* database management db.select: select values from table
r3.* 3D raster processing r3.stats: 3D raster statistics
t.* temporal data processing t.rast.aggregate: temporal aggregation
g.* general data management g.rename: renames map
d.* display d.rast: display raster map

These are the main groups of modules. There is few more for specific purposes. Note also that some modules have multiple dots in their names. This often suggests further grouping. For example, modules staring with v.net. deal with vector network analysis.

The name of the module helps to understand its function, for example v.in.lidar starts with v so it deals with vector maps, the name follows with in which indicates that the module is for importing the data into GRASS GIS Spatial Database and finally lidar indicates that it deals with lidar point clouds.

Finding and running a module

To find a module for your analysis, type the term into the search box into the Modules tab in the Layer Manager, then keep pressing Enter until you find your module.

Alternatively, you can just browse through the module tree in the Modules tab. You can also browse through the main menu. For example, to find information about a raster map, use: Raster → Reports and statistics → Basic raster metadata.

Running a module as a command

If you already know the name of the module, you can just use it in the command line. The GUI offers a Command console tab with command line specifically build for running GRASS GIS modules. If you type module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get GUI dialog for the module.

You can use the command line to run also whole commands for example when you get a command, i.e. module and list of parameters, in the instructions.

Command line vs. GUI interface

GRASS modules can be executed either through a GUI or command line interface. The GUI offers a user-friendly approach to executing modules where the user can navigate to data layers that they would like to analyze and modify processing options with simple check boxes. The GUI also offers an easily accessible manual on how to execute a model. The command line interface allows users to execute a module using command prompts specific to that module. This is handy when you are running similar analyses with minor modification or are familiar with the module commands for quick efficient processing. In this workshop we provide module prompts that can be copy and pasted into the command line for our workflow, but you can use both GUI and command line depending on personal preference. Look how GUI and command line interface represent the same tool.
Task: compute aspect (orientation) from provided digital elevation model using module r.slope.aspect using both module dialog and command line.

  • How to find modules? Modules are organized by their functionality in wxGUI menu, or we can search for them in Search modules tab. If we already know which module to use, we can just type it in the wxGUI command console.

Module parameters

The same analysis can be done using the following command:

r.neighbors -c input=elevation output=elev_smooth size=5

Conversely, you can fill the GUI dialog parameter by parameter when you have the command.

Computational region

Before we use a module to compute a new raster map, we must set properly computational region. All raster computations will be performed in the specified extent and with the given resolution.

Computational region is an important raster concept in GRASS GIS. In GRASS a computational region can be set, subsetting larger extent data for quicker testing of analysis or analysis of specific regions based on administrative units. We provide a few points to keep in mind when using the computational region function:

  • defined by region extent and raster resolution
  • applies to all raster operations
  • persists between GRASS sessions, can be different for different mapsets
  • advantages: keeps your results consistent, avoid clipping, for computationally demanding tasks set region to smaller extent, check your result is good and then set the computational region to the entire study area and rerun analysis
  • run g.region -p or in menu Settings - Region - Display region to see current region settings

The numeric values of computational region can be checked using:

g.region -p

After executing the command you will get something like this:

north:      220750
south:      220000
west:       638300
east:       639000
nsres:      1
ewres:      1
rows:       750
cols:       700
cells:      525000

Computational region can be set also using a vector map. In that case, only extent is set (as vector map does not have any resolution - at least not in the way raster map does). In GUI, this can be done in the same way as for the raster map. In the command line, it looks like this:

g.region vector=lakes

Resolution can be set separately using the res parameter of the g.region module. The units are the units of the current location, in our case meters. This can be done in the Resolution tab of the g.region dialog or in the command line in the following way (using also the -a flag to print the new values):

g.region res=3 -p

The new resolution may be slightly modified in this case to fit into the extent which we are not changing. However, often we want the resolution to be the exact value we provide and we are fine with a slight modification of the extent. That's what -a flag is for.

The following example command will use the extent from the vector named lakes, use resolution 10, modify the extent to align it to this 10 meter resolution, and print the values of this new computational region settings:

g.region vector=lakes res=10 -a -p

Running modules

Find the module for computing viewshed in the menu or the module tree under Raster → Terrain analysis → Visibility and aspect or simply run r.viewshed from the Console.

r.viewshed input=elevation output=viewshed coordinates=638587.6,220649.5

3D view

We can explore our study area in 3D view.

  1. Add elev_lid792_1m and uncheck or remove any other layers.
  2. Set computational region to this raster. Switch to 3D view (in the right corner on Map Display).
  3. Adjust the view (perspective, height, vertical exaggeration)
  4. In Data tab, set Fine mode resolution to 1 and set slope (computed in the previous step) as the color of the surface.
  5. When finished, switch back to 2D view.
3D visualization of elev_lid792_1m DEM with slope draped over

Introduction to scripting in GRASS GIS with Python

The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python editor integrated in GRASS GIS accessible from the toolbar or the Python tab in the Layer Manager. Another option is to write the Python code in your favorite plain text editor like Notepad++ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu File -> Launch script.



The GRASS GIS 7 Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include:

Besides, this library provides several wrapper functions for often called modules.

Calling GRASS GIS modules

We will use GRASS GUI Simple Python Editor to run the commands. You can open it from Python tab. For longer scripts, you can create a text file, save it into your current working directory and run it with python myscript.py from the GUI command console or terminal.

When you open Simple Python Editor, you find a short code snippet. It starts with importing GRASS GIS Python Scripting Library:

import grass.script as gscript

In the main function we call g.region to see the current computational region settings:

gscript.run_command('g.region', flags='p')

Note that the syntax is similar to bash syntax (g.region -p), just the flag is specified in a parameter. Now we can run the script by pressing the Run button in the toolbar. In Layer Manager we get the output of g.region.

Before running any GRASS raster modules, you need to set the computational region. In this example, we set the computational extent and resolution to the raster layer elevation. Replace the previous g.region command with the following line:

gscript.run_command('g.region', raster='elevation')

The run_command() function is the most commonly used one. We will use it to compute viewshed using r.viewshed. Add the following line after g.region command:

gscript.run_command('r.viewshed', input='elevation', output='viewshed_python', coordinates=TODO, overwrite=True)

Parameter overwrite is needed if we rerun the script and rewrite the raster. Now let's look at how big the viewshed is by using r.univar. Here we use parse_command to obtain the statistics as a Python dictionary

univar = gscript.parse_command('r.univar', map='viewshed_python', flags='g')
print univar['n']

The printed result is the number of cells of the viewshed.

List of convenient wrapper functions with examples:

Introduction to Blender

Animation in GRASS GIS

Scripting of rendering in GRASS GIS

import grass.script as gs

for point_id in [0.5, 1, 2, 5]:
    gs.run_command('d.mon', start='cairo',
        output='viewshed_{}.png'.format(point_id))
    gs.run_command('d.rast',
        map='viewshed_{}'.format(point_id))
    gs.run_command('d.mon', stop='cairo')

To

os.environ['GRASS_FONT'] = 'sans'
os.environ['GRASS_RENDER_IMMEDIATE'] = 'cairo'
os.environ['GRASS_RENDER_FILE_READ'] = 'TRUE'
os.environ['GRASS_LEGEND_FILE'] = 'legend.txt'
os.environ['GRASS_RENDER_FILE'] = 'map.png'

See also

GRASS GIS at ICC 2017 Unleash the power of GRASS GIS at US-IALE 2017