Analytical data visualizations at ICC 2017: Difference between revisions
⚠️Wenzeslaus (talk | contribs) |
mNo edit summary |
||
(97 intermediate revisions by 3 users not shown) | |||
Line 1: | Line 1: | ||
[[File:icc2017 logo.png|300px|right|none]] | [[File:icc2017 logo.png|300px|right|none]] | ||
[[File:Grassgis logo colorlogo text whitebg.png|300px|right|none]] | [[File:Grassgis logo colorlogo text whitebg.png|250px|right|none]] | ||
[[File:Blender-plain.png|300px|right|none]] | |||
This is material for two sessions at the [[GRASS GIS at ICC 2017|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 | This is material for two sessions at the [[GRASS GIS at ICC 2017|ICC 2017]] SDI-Open (Spatial data infrastructures, standards, open source and open data for geospatial) 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 and Blender, example of their processing capabilities and visualization techniques. Participants interactively visualize open data and design maps with several open source geospatial tools including Tangible Landscape system. | ||
Authors: [[User:Wenzeslaus|Vaclav Petras]], [[User:Annakrat|Anna Petrasova]], Payam Tabrizian, Brendan Harmon, and Helena Mitasova | |||
Testing and review: Garrett Millar | |||
We start the workshop with a [https://ncsu-geoforall-lab.github.io/ICC-workshop-presentation/index.html#/ presentation] about the overview of geospatial tools used in this workshop. | We start the workshop with a [https://ncsu-geoforall-lab.github.io/ICC-workshop-presentation/index.html#/ presentation] about the overview of geospatial tools used in this workshop. | ||
Note: Screenshots of GUI may be outdated, however the rest of the tutorial is updated to work with GRASS GIS 7.8. | |||
== Data == | == Data == | ||
=== Prepared data === | === Prepared data === | ||
[http://fatra.cnr.ncsu.edu/icc2017/ICC_workshop_data.zip Workshop dataset] contains: | |||
* digital surface model derived from North Carolina Flood Plain Mapping lidar | * digital surface model (1m resolution) derived from publicly available North Carolina Flood Plain Mapping lidar | ||
* orthophoto | * orthophoto (0.5m resolution) from USGS | ||
=== Download OpenStreetMap data === | === Download OpenStreetMap data === | ||
Line 16: | Line 24: | ||
[https://en.wikipedia.org/wiki/OpenStreetMap OpenStreetMap] data. | [https://en.wikipedia.org/wiki/OpenStreetMap OpenStreetMap] data. | ||
You can use Wizard to create simple queries on the currently zoomed extent. | You can use Wizard to create simple queries on the currently zoomed extent. For example, zoom to a small area and paste this into the Wizard and run the query: | ||
highway=* and type:way | highway=* and type:way | ||
The query was built in the editor and ran | The query was built in the editor and ran so you can now see the results in a map. | ||
Now, paste this query in the editor window: | Now, paste this query in the editor window: | ||
Line 27: | Line 35: | ||
// gather results | // gather results | ||
( | ( | ||
// query part for: | // query part for: "highway=*" | ||
way["highway"](35.76599,-78.66249,35.77230,-78.65261); | way["highway"](35.76599,-78.66249,35.77230,-78.65261); | ||
); | ); | ||
Line 44: | Line 52: | ||
// gather results | // gather results | ||
( | ( | ||
// query part for: | // query part for: "building=*" | ||
way["building"](35.76599,-78.66249,35.77260,-78.65261); | way["building"](35.76599,-78.66249,35.77260,-78.65261); | ||
relation["building"](35.76599,-78.66249,35.77260,-78.65261); | |||
); | ); | ||
// print results | // print results | ||
Line 56: | Line 64: | ||
In the next steps we will import this data into GRASS GIS. | In the next steps we will import this data into GRASS GIS. | ||
== Introduction to GRASS GIS == | == Session 1: Introduction to GRASS GIS == | ||
Here we provide an overview of [https://grass.osgeo.org 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. | Here we provide an overview of [https://grass.osgeo.org 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 ==== | ==== Structure of the GRASS GIS Spatial Database ==== | ||
GRASS uses unique database terminology and structure ([https://grass.osgeo.org/grass72/manuals/grass_database.html 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 | GRASS uses unique database terminology and structure ([https://grass.osgeo.org/grass72/manuals/grass_database.html 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 your 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. | * A '''GRASS GIS Spatial Database''' (''GRASS database'') consists of directory with specific Locations (projects) where data (data layers/maps) are stored. | ||
Line 76: | Line 84: | ||
Open Location Wizard with button ''New'' in the left part of the welcome screen. Select a name for the new location, select EPSG method and code 3358. | Open Location Wizard with button ''New'' in the left part of the welcome screen. Select a name for the new location, select EPSG method and code 3358. | ||
When the wizard is finished, the new location will be listed on the start-up screen. Select the new location and mapset PERMANENT and press ''Start GRASS session''. | When the wizard is finished, the new location will be listed on the start-up screen. Select the new location and mapset PERMANENT and press ''Start GRASS session''. | ||
<center><gallery perrow=3 widths=400 heights=270> | <center><gallery perrow=3 widths=400 heights=270> | ||
Line 86: | Line 92: | ||
Image:ICC_workshop_location_wizard_summary_7.2.png| Review summary page and confirm | Image:ICC_workshop_location_wizard_summary_7.2.png| Review summary page and confirm | ||
</gallery></center> | </gallery></center> | ||
==== Current working directory ==== | |||
A current working directory is a concept which is separate from the GRASS database, location and mapset discussed above. The current working directory is a directory (also called folder) where any program (no just GRASS GIS) writes and reads files unless a path to the file is provided. When using GRASS GIS GUI, the current working directory is usually not used. However, when using command line or doing scripting, the current working directory is important. | |||
Change your current working directory to some place where you have read and write permissions. For example, you can create a directory <tt>ICC</tt> in the <tt>Documents</tt>. The current working directory can be changed from the GUI using ''Settings → GRASS working environment → Change working directory'' or from the Console using the <tt>cd</tt> command. | |||
==== Importing data ==== | ==== Importing data ==== | ||
[[File:Import_raster_7.2.1.png|400px|thumb|right|Import raster data]] | [[File:Import_raster_7.2.1.png|400px|thumb|right|Import raster data]] | ||
In this step we import into GRASS GIS | In this step we import the provided data into GRASS GIS. In menu ''File - Import raster data'' select ''Common formats import'' and in the dialog browse to find <tt>dsm.tif</tt> and click button ''Import''. Repeat for raster file <tt>ground.tif</tt> and <tt>ortho.tif</tt>. | ||
Similarly, go to ''File - Import vector data - Common formats import'' and import files <tt>roads.geojson</tt> and <tt>buildings.geojson</tt>. '''Note that in this case we need to change the output name from default <tt>OGRGeoJSON</tt> to <tt>roads</tt> and <tt>buildings</tt>, respectively.''' This vector data has different CRS, so a dialog informing you about the need to reproject the data appears and we confirm it. | Similarly, go to ''File - Import vector data - Common formats import'' and import files <tt>roads.geojson</tt> and <tt>buildings.geojson</tt>. '''Note that in this case we need to change the output name from default <tt>OGRGeoJSON</tt> to <tt>roads</tt> and <tt>buildings</tt>, respectively.''' This vector data has different CRS, so a dialog informing you about the need to reproject the data appears and we confirm it. | ||
All the imported layers should be added to GUI, if not, add them manually. | All the imported layers should be added to GUI, if not, add them manually (see below). | ||
<br clear=all> | <br clear=all> | ||
==== Displaying and exploring data ==== | ==== Displaying and exploring data ==== | ||
The GUI interface allows you to display raster | The GUI interface allows you to display raster and 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. | ||
<center><gallery perrow=2 widths=500 heights=250> | <center><gallery perrow=2 widths=500 heights=250> | ||
Line 104: | Line 116: | ||
Image:Wxgui_add_legend_elev_lid.png|Add raster legend | Image:Wxgui_add_legend_elev_lid.png|Add raster legend | ||
Image:GRASS FUTURES wxGUI display.png|Layer Manager and Map Display overview. Annotations show how to add raster layer, query, add legend. | Image:GRASS FUTURES wxGUI display.png|Layer Manager and Map Display overview. Annotations show how to add raster layer, query, add legend. | ||
Image:wxGUI_metadata_rinfo.png|Show raster map metadata by right | Image:wxGUI_metadata_rinfo.png|Show raster map metadata by right clicking on layer | ||
</gallery></center> | </gallery></center> | ||
=== GRASS GIS modules === | === GRASS GIS modules overview === | ||
One of the advantages of GRASS is the diversity and number of modules that let you analyze all | One of the advantages of GRASS is the diversity and number of modules that let you analyze all manners of spatial and temporal data. GRASS GIS has over [https://grass.osgeo.org/grass72/manuals/full_index.html 500 different modules] in the core distribution and over [https://grass.osgeo.org/grass72/manuals/addons/ 300 addon modules] that can be used to prepare and analyze data. | ||
GRASS functionality is available through ''modules'' (tools, functions). Modules respect the following naming conventions: | GRASS functionality is available through ''modules'' (tools, functions). Modules respect the following naming conventions: | ||
Line 135: | Line 147: | ||
|}</center> | |}</center> | ||
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. | These are the main groups of modules. There is a 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. | 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. | 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''. | 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''. | ||
If you already know the name of the module, you can just use it in the command line. The GUI offers a Console tab with command line specifically built for running GRASS GIS modules. If you type the module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get GUI dialog for the module. | |||
<center><gallery perrow=3 widths=380 heights=190> | <center><gallery perrow=3 widths=380 heights=190> | ||
Image:wxGUI_module_search.png|Search for a module in module tree (searches in names, descriptions and keywords) | Image:wxGUI_module_search.png|Search for a module in module tree (searches in names, descriptions and keywords) | ||
Image:wxGUI_menu_rinfo.png|Modules can be also found in the main menu | Image:wxGUI_menu_rinfo.png|Modules can be also found in the main menu | ||
Image:wxGUI_console_completion.png|Automatic suggestions when typing name of the module: By typing prefix r. we make a list of modules starting with that prefix to show up. | Image:wxGUI_console_completion.png|Automatic suggestions when typing name of the module: By typing prefix r. we make a list of modules starting with that prefix to show up. | ||
</gallery></center> | </gallery></center> | ||
You can use the command line to run | You can also use the command line to run complete commands, i.e. module and list of parameters. These commands are what you usually find in the documentation or in this material. You can use them directly in the Console tab or you can use them as instructions how to fill the GUI dialog. | ||
=== 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 | [[File:Wxgui_module_parameters_r_neighbors.png|400px|thumb|right| Module dialog]] | ||
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 copied and pasted into the command line for our workflow, but you can use both the GUI and command line interface depending on personal preference. Look how {{cmd|wxGUI.modules|version=71|desc=GUI and command line interface}} represent the same tool. | |||
The same analysis can be done using the following command: | The same analysis can be done using the following command: | ||
Line 175: | Line 177: | ||
Conversely, you can fill the GUI dialog parameter by parameter when you have the command. | Conversely, you can fill the GUI dialog parameter by parameter when you have the command. | ||
<br clear=all> | |||
=== Computational region === | === Computational region === | ||
Before we use a module to compute a new raster map, we must set | Before we use a module to compute a new raster map, we must properly set the 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: | 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: | ||
Line 184: | Line 187: | ||
* applies to all raster operations | * applies to all raster operations | ||
* persists between GRASS sessions, can be different for different mapsets | * 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 | * advantages: keeps your results consistent, avoid clipping, for computationally demanding tasks set region to smaller extent, check that your result is good and then set the computational region to the entire study area and rerun analysis | ||
* run <code>g.region -p</code> or in menu ''Settings'' - ''Region'' - ''Display region'' to see current region settings | * run <code>g.region -p</code> or in menu ''Settings'' - ''Region'' - ''Display region'' to see current region settings | ||
<center><gallery perrow=3 widths=380 heights=190> | <center><gallery perrow=3 widths=380 heights=190> | ||
Image:Computational_region_two_rasters.png|Computational region concept: A raster with large extent (blue) is displayed as well as another raster with smaller extent (green). The computational region (red) is now set to match the smaller raster, so all the computations are limited to the smaller raster extent even if the input is the larger raster. (Not shown on the image: Also the resolution, not only the extent, matches the resolution of the smaller raster.) | Image:Computational_region_two_rasters.png|Computational region concept: A raster with large extent (blue) is displayed as well as another raster with smaller extent (green). The computational region (red) is now set to match the smaller raster, so all the computations are limited to the smaller raster extent even if the input is the larger raster. (Not shown on the image: Also the resolution, not only the extent, matches the resolution of the smaller raster.) | ||
Image:WxGUI set region.png|Simple ways to set computational region from GUI | Image:WxGUI set region.png|Simple ways to set computational region from GUI: On the left, set region to match raster map. On the right, select the highlighted option and then set region by drawing rectangle. | ||
Image:Wxgui_computational_region_set_from_raster.png|Set computational region (extent and resolution) to match a raster (''Layers'' tab in the ''Layer Manager'') | Image:Wxgui_computational_region_set_from_raster.png|Set computational region (extent and resolution) to match a raster (''Layers'' tab in the ''Layer Manager'') | ||
</gallery></center> | </gallery></center> | ||
Line 215: | Line 218: | ||
=== Running modules === | === Running modules === | ||
Find the module for computing viewshed in the menu or the module tree under ''Raster → Terrain analysis → Visibility'' | Find the module for computing viewshed in the menu or the module tree under ''Raster → Terrain analysis → Visibility'', or simply run {{cmd|r.viewshed}} from the ''Console''. | ||
r.viewshed input=dsm output= | r.viewshed input=dsm output=dsm_viewshed coordinates=640167,223907 observer_elevation=1.72 target_elevation=1 max_distance=400 | ||
=== 3D view === | === 3D view === | ||
We can explore our study area in 3D view | [[File:ICC_workshop_3Dview_ortho.png|800px|thumb|right| 3D visualization of DSM with orthophoto draped over]] | ||
# Add <tt> | We can explore our study area in 3D view (use image on the right if clarification is needed for below steps): | ||
# | # Add raster <tt>dsm</tt> and '''uncheck or remove''' any other layers. Any layer in Layer Manager is interpreted as surface in 3D view. | ||
# Adjust the view (perspective, height | # Switch to 3D view (in the right corner of Map Display). | ||
# In ''Data'' tab, set ''Fine mode resolution'' to 1 and set <tt> | # Adjust the view (perspective, height). | ||
# In ''Data'' tab, set ''Fine mode resolution'' to 1 and set <tt>ortho</tt> as the color of the surface (the orthophoto is draped over the DSM). | |||
# Go back to ''View'' tab and explore the different view directions using the green puck. | |||
# Go again to ''Data'' tab and change color to <tt>viewshed</tt> raster computed in the previous steps. | |||
# Go to ''Appearance'' tab and change the light conditions (lower the light height, change directions). | |||
# When finished, switch back to 2D view. | # When finished, switch back to 2D view. | ||
<br clear=all> | |||
== Introduction to scripting in GRASS GIS with Python == | == Session 1: 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 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. | ||
Line 240: | Line 247: | ||
The GRASS GIS 7 Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include: | 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}}''': | * '''{{pyapi|script|script.core|run_command}}''': 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|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|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 | * '''{{pyapi|script|script.core|write_command}}''': for modules expecting text input from either standard input or file | ||
We will use GRASS GUI Simple Python Editor to run the commands. You can open it from Python tab. | 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. | 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. | ||
Line 254: | Line 260: | ||
It starts with importing GRASS GIS Python Scripting Library: | It starts with importing GRASS GIS Python Scripting Library: | ||
<source lang="python"> | <source lang="python"> | ||
import grass.script as | import grass.script as gs | ||
</source> | </source> | ||
Line 260: | Line 266: | ||
<source lang="python"> | <source lang="python"> | ||
gs.run_command('g.region', flags='p') | |||
</source> | </source> | ||
Note that the syntax is similar to bash syntax (<code>g.region -p</code>), | Note that the syntax is similar to bash syntax (<code>g.region -p</code>), only 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>. | 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>. | ||
Line 269: | Line 275: | ||
<source lang="python"> | <source lang="python"> | ||
gs.run_command('g.region', raster='dsm') | |||
</source> | </source> | ||
Line 276: | Line 282: | ||
<source lang="python"> | <source lang="python"> | ||
gs.run_command('r.viewshed', input='dsm', output='python_viewshed', coordinates=(640645, 224035), overwrite=True) | |||
</source> | </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 | 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 | {{cmd|r.univar}}. Here we use <code>parse_command</code> to obtain the statistics as a Python dictionary | ||
<source lang="python"> | <source lang="python"> | ||
univar = | univar = gs.parse_command('r.univar', map='python_viewshed', flags='g') | ||
print univar['n'] | print(univar['n']) | ||
</source> | </source> | ||
The printed result is the number of cells of the viewshed. | The printed result is the number of cells of the viewshed. | ||
List of convenient wrapper functions with examples: | GRASS GIS Python Scripting Library also provides several wrapper functions for often called modules. List of convenient wrapper functions with examples includes: | ||
* Raster metadata using {{pyapi|script|script.raster|raster_info}}: <code> | * Raster metadata using {{pyapi|script|script.raster|raster_info}}: <code>gs.raster_info('viewshed_python')</code> | ||
* Vector metadata using {{pyapi|script|script.vector|vector_info}}: <code> | * Vector metadata using {{pyapi|script|script.vector|vector_info}}: <code>gs.vector_info('roads')</code> | ||
* List raster data in current location using {{pyapi|script|script.core|list_grouped}}: <code> | * List raster data in current location using {{pyapi|script|script.core|list_grouped}}: <code>gs.list_grouped(type=['raster'])</code> | ||
* Get current computational region using {{pyapi|script|script.core|region}}: <code> | * Get current computational region using {{pyapi|script|script.core|region}}: <code>gs.region()</code> | ||
Here are two commands (to be executed in the ''Console'') often used when working with scripts. First is setting the computational region. We can do that in a script, but it better and more general to do it before executing the script (so that the script can be used with different computational region settings): | Here are two commands (to be executed in the ''Console'') often used when working with scripts. First is setting the computational region. We can do that in a script, but it is better and more general to do it before executing the script (so that the script can be used with different computational region settings): | ||
g.region raster=dsm | g.region raster=dsm | ||
Line 305: | Line 313: | ||
The above command actually won't remove the maps, but it will inform you which it will remove if you provide the <tt>-f</tt> flag. | The above command actually won't remove the maps, but it will inform you which it will remove if you provide the <tt>-f</tt> flag. | ||
== Batch processing using Python == | == Session 1: Batch processing using Python == | ||
In this example, we will calculate the viewsheds along a road to simulate in a simple way what a driver would see. | In this example, we will calculate the viewsheds along a road to simulate, in a simple way, what a driver would see. First we will prepare view points by running the modules directly through the GUI dialogs or in the Console. Then we will use Python to process all the created points. | ||
First we extract only segments of road with name 'Umstead Drive': | First we extract only segments of road with name 'Umstead Drive' (use the code in the Console tab or use the GUI dialog): | ||
v.extract input=roads where="name = 'Umstead Drive'" output=umstead_drive_segments | v.extract input=roads where="name = 'Umstead Drive'" output=umstead_drive_segments | ||
We will join the segments into one polyline using {{cmd|v.build.polylines}} | We will join the segments into one polyline using {{cmd|v.build.polylines}} | ||
v.build.polylines input=umstead_drive_segments output=umstead_drive cats=first | v.build.polylines input=umstead_drive_segments output=umstead_drive cats=first | ||
Now we generate viewpoints in regular intervals along the line: | Now we generate viewpoints in regular intervals along the line: | ||
v.to.points input=umstead_drive type=line output=viewpoints dmax=50 | v.to.points input=umstead_drive type=line output=viewpoints dmax=50 | ||
We compute viewshed from each point using a Python script. | We compute viewshed from each point using a Python script. | ||
We will use the difference between the DSM and ground elevation to filter out points that happen to lie on some vegetation covering the road. This is to avoid computing viewsheds from the top of a tree, | |||
which would significantly distort viewshed size. | |||
r.mapcalc "diff = dsm - ground" | |||
r.colors map=diff color=differences | |||
From these viewsheds we can then easily compute ''Cumulative viewshed'' (how much a place is visible). | From these viewsheds we can then easily compute ''Cumulative viewshed'' (how much a place is visible). | ||
For that we will use module {{cmd|r.series}}, which overlays the viewshed rasters and computes how many times each cell is visible in all the rasters. | For that we will use module {{cmd|r.series}}, which overlays the viewshed rasters and computes how many times each cell is visible in all the rasters. | ||
Copy the following code | Copy the following code | ||
into Python Simple Editor and press Run button. | into Python Simple Editor (be sure to overwrite any already existing code) and press Run button. | ||
<source lang="python"> | <source lang="python"> | ||
#!/usr/bin/env | #!/usr/bin/env python3 | ||
import grass.script as | import grass.script as gs | ||
def main(): | def main(): | ||
# obtain the coordinates of viewpoints | # obtain the coordinates of viewpoints | ||
viewpoints = | viewpoints = gs.read_command('v.out.ascii', input='viewpoints', | ||
separator='comma', layer=2).strip() | |||
# loop through the viewpoints and compute viewshed from each of them | # loop through the viewpoints and compute viewshed from each of them | ||
for point in viewpoints.splitlines(): | for point in viewpoints.splitlines(): | ||
if not point: | |||
# skip empty lines | |||
continue | |||
x, y, cat = point.split(',') | x, y, cat = point.split(',') | ||
diff = gs.raster_what(map='diff', coord=(float(x), float(y)))[0]['diff']['value'] | |||
if float(diff) > 2: | |||
continue | |||
gs.run_command('r.viewshed', input='dsm', output='viewshed' + cat, | |||
coordinates=(x, y), max_distance=300, overwrite=True) | |||
# obtain all viewshed results and set their color to yellow | # obtain all viewshed results and set their color to yellow | ||
# export the list of viewshed map names to a file | # export the list of viewshed map names to a file | ||
maps_file = | maps_file = 'viewsheds.txt' | ||
gs.run_command('g.list', type='raster', pattern='viewshed*', output=maps_file) | |||
gs.write_command('r.colors', file=maps_file, rules='-', | |||
stdin='0% yellow \n 100% yellow') | |||
# cumulative viewshed | # cumulative viewshed | ||
gs.run_command('r.series', file='viewsheds.txt', output='cumulative_viewshed', method='count') | |||
# set color of cumulative viewshed from grey to yellow | # set color of cumulative viewshed from grey to yellow | ||
gs.write_command('r.colors', map='cumulative_viewshed', rules='-', | |||
stdin='0% 70:70:70 \n 100% yellow') | |||
if __name__ == '__main__': | if __name__ == '__main__': | ||
Line 359: | Line 377: | ||
</source> | </source> | ||
We can visualize the cumulative viewshed draped over the DSM using 3D view. | We can visualize the cumulative viewshed draped over the DSM using 3D view. If needed review the instructions in [[#3D_view| previous section about 3D view]]. | ||
[[File:cumulative_viewshed.png|800px|thumb|center| Cumulative viewshed computed along a road]] | [[File:cumulative_viewshed.png|800px|thumb|center| Cumulative viewshed computed along a road]] | ||
To proceed with the next part, we will export the viewpoints and cumulative viewshed in order for Blender to use it. | To proceed with the next part, we will export the viewpoints and cumulative viewshed in order for Blender to use it. | ||
Copy and paste the individual lines to GRASS GUI command console and execute them separately, omit the lines starting with #. | Copy and paste the individual lines to GRASS GUI command console and execute them '''separately''', omit the lines starting with #. | ||
<source> | <source lang="bash"> | ||
# turn 2D points into 3D using the DSM information | # turn 2D points into 3D using the DSM information | ||
v.drape input=viewpoints output=viewpoints_3d elevation=dsm method=bilinear | v.drape input=viewpoints output=viewpoints_3d elevation=dsm method=bilinear | ||
Line 375: | Line 394: | ||
</source> | </source> | ||
== Introduction to Blender == | == Session 1: Introduction to Blender == | ||
Materials for the Blender part of this workshop are available here: | |||
== Animation in GRASS GIS == | https://github.com/ptabriz/ICC_2017_Workshop/blob/master/README.md | ||
== Session 2: Animation in GRASS GIS == | |||
In the following exercise we will visualize the viewsheds along the road (we computed above) as an animation. We will use {{cmd|g.gui.animation|desc=GRASS GIS Animation Tool}}. Start it from menu ''File - Animation Tool'': | In the following exercise we will visualize the viewsheds along the road (we computed above) as an animation. We will use {{cmd|g.gui.animation|desc=GRASS GIS Animation Tool}}. Start it from menu ''File - Animation Tool'': | ||
[[File:futures_anim_start_animation_tool.png|thumb|none|Animation tool in the main menu]] | [[File:futures_anim_start_animation_tool.png|thumb|none|Animation tool in the main menu]] | ||
Line 383: | Line 405: | ||
Start with ''Add new animation'' and click on ''Add space-time dataset or series of map layers''. In the ''Add space-time dataset layer'' click on a button next to the entry field and type <tt>^viewshed</tt>. | Start with ''Add new animation'' and click on ''Add space-time dataset or series of map layers''. In the ''Add space-time dataset layer'' click on a button next to the entry field and type <tt>^viewshed</tt>. | ||
This filters the rasters which start with "viewshed". | This filters the rasters which start with "viewshed". You can use [https://docs.python.org/2/library/re.html Python regular expressions] to filter maps in this dialog. | ||
Check only the desired raster maps are checked and then confirm the dialog and return to the dialog ''Add new animation''. | |||
Next we want to add orthophoto as base layer. Use ''Add raster map layer'' and select raster ortho. | Next we want to add orthophoto as base layer. Use ''Add raster map layer'' and select raster <tt>ortho</tt>. | ||
{|style="margin: 0 auto;" | {|style="margin: 0 auto;" | ||
Line 393: | Line 416: | ||
Do the same for the road vector <tt>umstead_drive</tt>, just instead of raster, add vector. Then rearrange the layers so that the ortho is below all other layers. | Do the same for the road vector <tt>umstead_drive</tt>, just instead of raster, add vector. Then, rearrange the layers so that the ortho is below all the other layers. | ||
Confirm the dialog and the animation should be loaded. | Confirm the dialog and the animation should be then loaded. | ||
{|style="margin: 0 auto;" | {|style="margin: 0 auto;" | ||
| [[File:Animation_tool_layers.png|thumb|none|500px|Change the order of layers]] | | [[File:Animation_tool_layers.png|thumb|none|500px|Change the order of layers]] | ||
Line 403: | Line 426: | ||
* If you want to zoom in, go back to Map Display, and set the computational region to the extent you want using toolbar button 'Various zoom options - Set computational region extent interactively'. Then press ''Render map'' in the Animation tool toolbar and it will re-render the animation based on the current computational region extent. | * If you want to zoom in, go back to Map Display, and set the computational region to the extent you want using toolbar button 'Various zoom options - Set computational region extent interactively'. Then press ''Render map'' in the Animation tool toolbar and it will re-render the animation based on the current computational region extent. | ||
== Scripting of rendering in GRASS GIS == | == Session 2: Scripting of rendering in GRASS GIS == | ||
Rendering of maps can be done in GRASS GIS using Python API in way which follows what is done in the GUI. Similarly to running the processing modules, GUI dialogs have a ''Copy'' button which will put a command line equivalent of what is | Rendering of maps can be done in GRASS GIS using Python API in a way which follows what is done in the GUI. Similarly to running the processing modules, GUI dialogs have a ''Copy'' button which will put a command line equivalent of what is currently set in the GUI dialog. For example, a raster properties dialog can give the code: <code>d.rast map=viewshed</code>. This code is written in Python syntax as <code>run_command('d.rast', map='viewshed')</code>. In addition to the rendering commands, we need to specify where (filename) and how (e.g. file format) to render this, which can be done using the {{cmd|d.mon}} command as below or environmental variables (in GUI we don't need to say anything because the context is clear). A full example is shown below: | ||
<source lang="python"> | <source lang="python"> | ||
import grass.script as | #!/usr/bin/env python | ||
import grass.script as gs | |||
viewsheds = | viewsheds = gs.list_strings(type=['raster'], pattern='viewshed*') | ||
for viewshed in viewsheds: | for viewshed in viewsheds: | ||
# drop mapset name from the image file name | |||
filename = '{}.png'.format(viewshed.split('@')[0]) | |||
# begin rendering | |||
gs.run_command('d.mon', start='cairo', output=filename) | |||
# rendering commands | |||
gs.run_command('d.rast', map='ortho') | |||
gs.run_command('d.vect', map='umstead_drive', color='232:232:232', width=3) | |||
gs.run_command('d.rast', map=viewshed) | |||
# end rendering | |||
gs.run_command('d.mon', stop='cairo') | |||
</source> | </source> | ||
In scripting, the {{cmd|d.mon}} is actually often not used because only one monitor can be active. This is very convenient simplification for interactive use, but limiting for scripting. When the rendering runs in parallel, it is necessary to use environmental {{cmd|variables}} instead of the {{cmd|d.mon}} module. Variables are set globally (using <code>os.environ</code> dictionary) or for each function call (using <code>env</code> parameter). The following three variables are the most important: <tt>GRASS_RENDER_IMMEDIATE</tt> (usually set to <tt>cairo</tt>), <tt>GRASS_RENDER_FILE_READ</tt> (set to <tt>TRUE</tt>), and <tt>GRASS_RENDER_FILE</tt> (resulting image filename). Additional variables include <tt>GRASS_FONT</tt> (name of the font to be used, e.g. sans) and <tt>GRASS_LEGEND_FILE</tt> (file for vector legend, can be path to a temporary file). | The numbers <code>232:232:232</code> in the {{cmd|d.vect}} call are a GRASS RGB triplet. In other words, red, green, and blue colors are expressed as numbers from 0 to 255 separated by a colon (using HTML syntax, an equivalent would be <code>rgb(232,232,232)</code> or <code>#E8E8E8</code>). | ||
'''Information for scripting, parallelization, and supercomputing:''' In scripting, the {{cmd|d.mon}} is actually often not used because only one monitor can be active. This is a very convenient simplification for interactive use, but limiting for scripting. When the rendering runs in parallel, it is necessary to use environmental {{cmd|variables}} instead of the {{cmd|d.mon}} module. Variables are set globally (using <code>os.environ</code> dictionary) or for each function call (using <code>env</code> parameter). The following three variables are the most important: <tt>GRASS_RENDER_IMMEDIATE</tt> (usually set to <tt>cairo</tt>), <tt>GRASS_RENDER_FILE_READ</tt> (set to <tt>TRUE</tt>), and <tt>GRASS_RENDER_FILE</tt> (resulting image filename). Additional variables include: <tt>GRASS_FONT</tt> (name of the font to be used, e.g. sans) and <tt>GRASS_LEGEND_FILE</tt> (file for vector legend, can be path to a temporary file). | |||
'''Further exploration:''' | |||
* Use the <tt>roads</tt> vector map instead of just <tt>umstead_drive</tt>. | |||
* Add the <tt>buildings</tt> vector map in the same way as the <tt>roads</tt>. This time, you need to set two color options: color (for outline) and fill_color (for inner part). Use the vector map layer properties dialog in GUI to create the desired styling. Use the Copy button to get parameters in the command line syntax and rewrite it to Python. You can also just use the following colors: 188:82:47, 213:115:82 (#BC522F, #D57352) or 170:166:157, 213:205:189 (#AAA69D, #D5CDBD). | |||
* If you have ImageMagic available, create an animated GIF using something like: <code>convert 'viewshed?.png' 'viewshed??.png' viewshed.gif</code> | |||
== Session 2: Light up the terrain with viewsheds == | |||
Materials for the Blender part of this workshop are available here: | |||
https://github.com/ptabriz/ICC_2017_Workshop/blob/master/README.md#light-up-the-terrain-with-viewsheds | |||
== Session 2: Tangible Landscape == | |||
''Tangibly explore landscapes'' with [http://tangible-landscape.github.io/ Tangible Landscape] | |||
[[File:Ncgis2017 ncsu booth tangible landscape.jpg|300px|thumb|right|Tangible Landscape (photos by Anna Klevtcova)]] | |||
== Software == | == Software == | ||
Line 427: | Line 477: | ||
'''MS Windows''' | '''MS Windows''' | ||
Download the standalone GRASS GIS binaries ([https://grass.osgeo.org/grass72/binary/mswindows/native/x86_64/WinGRASS-7.2.1-1-Setup-x86_64.exe 64-bit version], or [https://grass.osgeo.org/grass72/binary/mswindows/native/x86/WinGRASS-7.2.1-1-Setup-x86.exe 32-bit version]) from [https://grass.osgeo.org/ grass.osgeo.org]. | |||
'''Mac OS''' | '''Mac OS''' | ||
Install GRASS GIS | Install GRASS GIS using [https://brew.sh/ Homebrew] [https://github.com/OSGeo/homebrew-osgeo4mac osgeo4mac]: | ||
<pre> | <pre> | ||
brew tap osgeo/osgeo4mac | brew tap osgeo/osgeo4mac | ||
Line 452: | Line 501: | ||
''' Linux ''' | ''' Linux ''' | ||
For other Linux distributions other then Ubuntu, please try to find GRASS GIS in their package managers. | For other Linux distributions other then Ubuntu, please try to find GRASS GIS in their package managers. | ||
Line 469: | Line 511: | ||
* [[Introduction to GRASS GIS with terrain analysis examples]] | * [[Introduction to GRASS GIS with terrain analysis examples]] | ||
* [[GRASS Location Wizard]] | * [[GRASS Location Wizard]] | ||
* [[GRASS and Blender]] | |||
[[Category: Tutorial]] | [[Category: Tutorial]] | ||
[[Category: Workshops]] | [[Category: Workshops]] | ||
[[Category: 2017]] | [[Category: 2017]] |
Latest revision as of 18:10, 14 September 2020
This is material for two sessions at the ICC 2017 SDI-Open (Spatial data infrastructures, standards, open source and open data for geospatial) 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 and Blender, example of their processing capabilities and visualization techniques. Participants interactively visualize open data and design maps with several open source geospatial tools including Tangible Landscape system.
Authors: Vaclav Petras, Anna Petrasova, Payam Tabrizian, Brendan Harmon, and Helena Mitasova
Testing and review: Garrett Millar
We start the workshop with a presentation about the overview of geospatial tools used in this workshop.
Note: Screenshots of GUI may be outdated, however the rest of the tutorial is updated to work with GRASS GIS 7.8.
Data
Prepared data
Workshop dataset contains:
- digital surface model (1m resolution) derived from publicly available North Carolina Flood Plain Mapping lidar
- orthophoto (0.5m resolution) from USGS
Download OpenStreetMap data
We will use overpass turbo (web-based data filtering tool) to create and run Overpass API queries to obtain OpenStreetMap data.
You can use Wizard to create simple queries on the currently zoomed extent. For example, zoom to a small area and paste this into the Wizard and run the query:
highway=* and type:way
The query was built in the editor and ran so you can now see the results in a map.
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). Name the vector file roads.geojson.
Now let's also download all building footprints and again export them as GeoJSON. Name the vector file buildings.geojson.
[out:json][timeout:25]; // gather results ( // query part for: "building=*" way["building"](35.76599,-78.66249,35.77260,-78.65261); relation["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.
Session 1: 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 your 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.
-
GRASS GIS Spatial Database structure
Creating a GRASS database for the tutorial
Start GRASS GIS, a start-up screen should appear. 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.
We will create a new location for our project with CRS (coordinate reference system) NC State Plane Meters with EPSG code 3358. Open Location Wizard with button New in the left part of the welcome screen. Select a name for the new location, select EPSG method and code 3358. When the wizard is finished, the new location will be listed on the start-up screen. Select the new location and mapset PERMANENT and press Start GRASS session.
-
GRASS GIS 7.2 startup dialog
-
Start Location Wizard and type the new location's name
-
Select method for describing CRS
-
Find and select EPSG 3358
-
Review summary page and confirm
Current working directory
A current working directory is a concept which is separate from the GRASS database, location and mapset discussed above. The current working directory is a directory (also called folder) where any program (no just GRASS GIS) writes and reads files unless a path to the file is provided. When using GRASS GIS GUI, the current working directory is usually not used. However, when using command line or doing scripting, the current working directory is important.
Change your current working directory to some place where you have read and write permissions. For example, you can create a directory ICC in the Documents. The current working directory can be changed from the GUI using Settings → GRASS working environment → Change working directory or from the Console using the cd command.
Importing data
In this step we import the provided data into GRASS GIS. In menu File - Import raster data select Common formats import and in the dialog browse to find dsm.tif and click button Import. Repeat for raster file ground.tif and ortho.tif.
Similarly, go to File - Import vector data - Common formats import and import files roads.geojson and buildings.geojson. Note that in this case we need to change the output name from default OGRGeoJSON to roads and buildings, respectively. This vector data has different CRS, so a dialog informing you about the need to reproject the data appears and we confirm it.
All the imported layers should be added to GUI, if not, add them manually (see below).
Displaying and exploring data
The GUI interface allows you to display raster and 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.
-
Add raster map layer
-
Add raster legend
-
Layer Manager and Map Display overview. Annotations show how to add raster layer, query, add legend.
-
Show raster map metadata by right clicking on layer
GRASS GIS modules overview
One of the advantages of GRASS is the diversity and number of modules that let you analyze all manners of spatial and temporal data. GRASS GIS has over 500 different modules in the core distribution and over 300 addon modules that can be used to prepare and analyze data.
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 a 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.
If you already know the name of the module, you can just use it in the command line. The GUI offers a Console tab with command line specifically built for running GRASS GIS modules. If you type the module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get GUI dialog for the module.
-
Search for a module in module tree (searches in names, descriptions and keywords)
-
Modules can be also found in the main menu
-
Automatic suggestions when typing name of the module: By typing prefix r. we make a list of modules starting with that prefix to show up.
You can also use the command line to run complete commands, i.e. module and list of parameters. These commands are what you usually find in the documentation or in this material. You can use them directly in the Console tab or you can use them as instructions how to fill the GUI dialog.
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 copied and pasted into the command line for our workflow, but you can use both the GUI and command line interface depending on personal preference. Look how GUI and command line interface represent the same tool.
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 properly set the 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 that 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
-
Computational region concept: A raster with large extent (blue) is displayed as well as another raster with smaller extent (green). The computational region (red) is now set to match the smaller raster, so all the computations are limited to the smaller raster extent even if the input is the larger raster. (Not shown on the image: Also the resolution, not only the extent, matches the resolution of the smaller raster.)
-
Simple ways to set computational region from GUI: On the left, set region to match raster map. On the right, select the highlighted option and then set region by drawing rectangle.
-
Set computational region (extent and resolution) to match a raster (Layers tab in the Layer Manager)
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 using a raster map:
g.region raster=dsm -p
Running modules
Find the module for computing viewshed in the menu or the module tree under Raster → Terrain analysis → Visibility, or simply run r.viewshed from the Console.
r.viewshed input=dsm output=dsm_viewshed coordinates=640167,223907 observer_elevation=1.72 target_elevation=1 max_distance=400
3D view
We can explore our study area in 3D view (use image on the right if clarification is needed for below steps):
- Add raster dsm and uncheck or remove any other layers. Any layer in Layer Manager is interpreted as surface in 3D view.
- Switch to 3D view (in the right corner of Map Display).
- Adjust the view (perspective, height).
- In Data tab, set Fine mode resolution to 1 and set ortho as the color of the surface (the orthophoto is draped over the DSM).
- Go back to View tab and explore the different view directions using the green puck.
- Go again to Data tab and change color to viewshed raster computed in the previous steps.
- Go to Appearance tab and change the light conditions (lower the light height, change directions).
- When finished, switch back to 2D view.
Session 1: 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.
-
Simple Python Editor integrated in GRASS GIS (since version 7.2) with Python tab in the background which contains an interactive Python shell.
-
Python tab with an interactive Python shell
The GRASS GIS 7 Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include:
- script.core.run_command(): used with modules which output raster/vector data where text output is not expected
- script.core.read_command(): used when we are interested in text output
- script.core.parse_command(): used with modules producing text output as key=value pair
- script.core.write_command(): for modules expecting text input from either standard input or file
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 gs
In the main function we call g.region to see the current computational region settings:
gs.run_command('g.region', flags='p')
Note that the syntax is similar to bash syntax (g.region -p
), only 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:
gs.run_command('g.region', raster='dsm')
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:
gs.run_command('r.viewshed', input='dsm', output='python_viewshed', coordinates=(640645, 224035), 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 = gs.parse_command('r.univar', map='python_viewshed', flags='g')
print(univar['n'])
The printed result is the number of cells of the viewshed.
GRASS GIS Python Scripting Library also provides several wrapper functions for often called modules. List of convenient wrapper functions with examples includes:
- Raster metadata using script.raster.raster_info():
gs.raster_info('viewshed_python')
- Vector metadata using script.vector.vector_info():
gs.vector_info('roads')
- List raster data in current location using script.core.list_grouped():
gs.list_grouped(type=['raster'])
- Get current computational region using script.core.region():
gs.region()
Here are two commands (to be executed in the Console) often used when working with scripts. First is setting the computational region. We can do that in a script, but it is better and more general to do it before executing the script (so that the script can be used with different computational region settings):
g.region raster=dsm
The second command is handy when we want to run the script again and again. In that case, we first need to remove the created raster maps, for example:
g.remove type=raster pattern="viewshed*"
The above command actually won't remove the maps, but it will inform you which it will remove if you provide the -f flag.
Session 1: Batch processing using Python
In this example, we will calculate the viewsheds along a road to simulate, in a simple way, what a driver would see. First we will prepare view points by running the modules directly through the GUI dialogs or in the Console. Then we will use Python to process all the created points.
First we extract only segments of road with name 'Umstead Drive' (use the code in the Console tab or use the GUI dialog):
v.extract input=roads where="name = 'Umstead Drive'" output=umstead_drive_segments
We will join the segments into one polyline using v.build.polylines
v.build.polylines input=umstead_drive_segments output=umstead_drive cats=first
Now we generate viewpoints in regular intervals along the line:
v.to.points input=umstead_drive type=line output=viewpoints dmax=50
We compute viewshed from each point using a Python script. We will use the difference between the DSM and ground elevation to filter out points that happen to lie on some vegetation covering the road. This is to avoid computing viewsheds from the top of a tree, which would significantly distort viewshed size.
r.mapcalc "diff = dsm - ground" r.colors map=diff color=differences
From these viewsheds we can then easily compute Cumulative viewshed (how much a place is visible). For that we will use module r.series, which overlays the viewshed rasters and computes how many times each cell is visible in all the rasters.
Copy the following code into Python Simple Editor (be sure to overwrite any already existing code) and press Run button.
#!/usr/bin/env python3
import grass.script as gs
def main():
# obtain the coordinates of viewpoints
viewpoints = gs.read_command('v.out.ascii', input='viewpoints',
separator='comma', layer=2).strip()
# loop through the viewpoints and compute viewshed from each of them
for point in viewpoints.splitlines():
if not point:
# skip empty lines
continue
x, y, cat = point.split(',')
diff = gs.raster_what(map='diff', coord=(float(x), float(y)))[0]['diff']['value']
if float(diff) > 2:
continue
gs.run_command('r.viewshed', input='dsm', output='viewshed' + cat,
coordinates=(x, y), max_distance=300, overwrite=True)
# obtain all viewshed results and set their color to yellow
# export the list of viewshed map names to a file
maps_file = 'viewsheds.txt'
gs.run_command('g.list', type='raster', pattern='viewshed*', output=maps_file)
gs.write_command('r.colors', file=maps_file, rules='-',
stdin='0% yellow \n 100% yellow')
# cumulative viewshed
gs.run_command('r.series', file='viewsheds.txt', output='cumulative_viewshed', method='count')
# set color of cumulative viewshed from grey to yellow
gs.write_command('r.colors', map='cumulative_viewshed', rules='-',
stdin='0% 70:70:70 \n 100% yellow')
if __name__ == '__main__':
main()
We can visualize the cumulative viewshed draped over the DSM using 3D view. If needed review the instructions in previous section about 3D view.
To proceed with the next part, we will export the viewpoints and cumulative viewshed in order for Blender to use it. Copy and paste the individual lines to GRASS GUI command console and execute them separately, omit the lines starting with #.
# turn 2D points into 3D using the DSM information
v.drape input=viewpoints output=viewpoints_3d elevation=dsm method=bilinear
# export them, ignore errors coming from too long attribute column names
v.out.ogr input=viewpoints_3d output=viewpoints_3d.shp format=ESRI_Shapefile
# export cumulative viewshed as PNG
r.out.gdal input=cumulative_viewshed output=cumulative_viewshed.png format=PNG
Session 1: Introduction to Blender
Materials for the Blender part of this workshop are available here:
https://github.com/ptabriz/ICC_2017_Workshop/blob/master/README.md
Session 2: Animation in GRASS GIS
In the following exercise we will visualize the viewsheds along the road (we computed above) as an animation. We will use GRASS GIS Animation Tool. Start it from menu File - Animation Tool:
Start with Add new animation and click on Add space-time dataset or series of map layers. In the Add space-time dataset layer click on a button next to the entry field and type ^viewshed.
This filters the rasters which start with "viewshed". You can use Python regular expressions to filter maps in this dialog.
Check only the desired raster maps are checked and then confirm the dialog and return to the dialog Add new animation.
Next we want to add orthophoto as base layer. Use Add raster map layer and select raster ortho.
Do the same for the road vector umstead_drive, just instead of raster, add vector. Then, rearrange the layers so that the ortho is below all the other layers.
Confirm the dialog and the animation should be then loaded.
Further exploration:
- If you want to zoom in, go back to Map Display, and set the computational region to the extent you want using toolbar button 'Various zoom options - Set computational region extent interactively'. Then press Render map in the Animation tool toolbar and it will re-render the animation based on the current computational region extent.
Session 2: Scripting of rendering in GRASS GIS
Rendering of maps can be done in GRASS GIS using Python API in a way which follows what is done in the GUI. Similarly to running the processing modules, GUI dialogs have a Copy button which will put a command line equivalent of what is currently set in the GUI dialog. For example, a raster properties dialog can give the code: d.rast map=viewshed
. This code is written in Python syntax as run_command('d.rast', map='viewshed')
. In addition to the rendering commands, we need to specify where (filename) and how (e.g. file format) to render this, which can be done using the d.mon command as below or environmental variables (in GUI we don't need to say anything because the context is clear). A full example is shown below:
#!/usr/bin/env python
import grass.script as gs
viewsheds = gs.list_strings(type=['raster'], pattern='viewshed*')
for viewshed in viewsheds:
# drop mapset name from the image file name
filename = '{}.png'.format(viewshed.split('@')[0])
# begin rendering
gs.run_command('d.mon', start='cairo', output=filename)
# rendering commands
gs.run_command('d.rast', map='ortho')
gs.run_command('d.vect', map='umstead_drive', color='232:232:232', width=3)
gs.run_command('d.rast', map=viewshed)
# end rendering
gs.run_command('d.mon', stop='cairo')
The numbers 232:232:232
in the d.vect call are a GRASS RGB triplet. In other words, red, green, and blue colors are expressed as numbers from 0 to 255 separated by a colon (using HTML syntax, an equivalent would be rgb(232,232,232)
or #E8E8E8
).
Information for scripting, parallelization, and supercomputing: In scripting, the d.mon is actually often not used because only one monitor can be active. This is a very convenient simplification for interactive use, but limiting for scripting. When the rendering runs in parallel, it is necessary to use environmental variables instead of the d.mon module. Variables are set globally (using os.environ
dictionary) or for each function call (using env
parameter). The following three variables are the most important: GRASS_RENDER_IMMEDIATE (usually set to cairo), GRASS_RENDER_FILE_READ (set to TRUE), and GRASS_RENDER_FILE (resulting image filename). Additional variables include: GRASS_FONT (name of the font to be used, e.g. sans) and GRASS_LEGEND_FILE (file for vector legend, can be path to a temporary file).
Further exploration:
- Use the roads vector map instead of just umstead_drive.
- Add the buildings vector map in the same way as the roads. This time, you need to set two color options: color (for outline) and fill_color (for inner part). Use the vector map layer properties dialog in GUI to create the desired styling. Use the Copy button to get parameters in the command line syntax and rewrite it to Python. You can also just use the following colors: 188:82:47, 213:115:82 (#BC522F, #D57352) or 170:166:157, 213:205:189 (#AAA69D, #D5CDBD).
- If you have ImageMagic available, create an animated GIF using something like:
convert 'viewshed?.png' 'viewshed??.png' viewshed.gif
Session 2: Light up the terrain with viewsheds
Materials for the Blender part of this workshop are available here:
Session 2: Tangible Landscape
Tangibly explore landscapes with Tangible Landscape
Software
Software should be pre-installed at ICC, but the following instructions can be used to install software on participants' laptops.
We use GRASS GIS 7.2 and Blender 2.78.
MS Windows
Download the standalone GRASS GIS binaries (64-bit version, or 32-bit version) from grass.osgeo.org.
Mac OS
Install GRASS GIS 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
For other Linux distributions other then Ubuntu, please try to find GRASS GIS in their package managers.