Skip to content

Instantly share code, notes, and snippets.

@joergsteinkamp

joergsteinkamp/Grass-GIS.ipynb

Last active December 14, 2021 07:40
Show Gist options
  • Star 4 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save joergsteinkamp/43fb13a82aeab9b5c0dcc4433f97f1a6 to your computer and use it in GitHub Desktop.
Save joergsteinkamp/43fb13a82aeab9b5c0dcc4433f97f1a6 to your computer and use it in GitHub Desktop.
Download ZIP
Simple Opensource GIS examples with Grass and R
Display the source blob
Display the rendered blob
Raw
Grass-GIS.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Importing Arc/Info raster and MS Access DB into Grass-GIS\n",
"\n",
"The [ISRIC-WISE global soil database](http://data.isric.org/geonetwork/srv/eng/catalog.search#/metadata/d9eca770-29a4-4d95-bf93-f32e1ab419c3) consists of\n",
"- Raster basemap of unique IDs\n",
"- Database of various soil properties\n",
"\n",
"Here: soil depth \n",
"\n",
"## Prerequisites\n",
"\n",
"1. [Grass-GIS](https://grass.osgeo.org/)\n",
"2. [mdbtools](https://sourceforge.net/projects/mdbtools/)\n",
"3. [sqlite3](https://www.sqlite.org/)\n",
"\n",
"## Setting up the Grass environment \n",
"\n",
"- EPSG code 4326: geographic lat/lon projection (ellipsoid WGS84)\n",
"- global terrestrial extent at 0.5° resolution (without Antarctica)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[H\u001b[2Jprojection: 3 (Latitude-Longitude)\n",
"zone: 0\n",
"datum: wgs84\n",
"ellipsoid: wgs84\n",
"north: 83N\n",
"south: 56S\n",
"west: 180W\n",
"east: 180E\n",
"nsres: 0:30\n",
"ewres: 0:30\n",
"rows: 278\n",
"cols: 720\n",
"cells: 200160\n",
"\u001b[H\u001b[2J"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cleaning up temporary files...\n",
"access: No such file or directory\n",
"ERROR: LOCATION </Users/jsteinkamp/Grass/wise> not available\n",
"Creating new GRASS GIS location/mapset...\n",
"\n",
" __________ ___ __________ _______________\n",
" / ____/ __ \\/ | / ___/ ___/ / ____/ _/ ___/\n",
" / / __/ /_/ / /| | \\__ \\\\_ \\ / / __ / / \\__ \\\n",
" / /_/ / _, _/ ___ |___/ /__/ / / /_/ // / ___/ /\n",
" \\____/_/ |_/_/ |_/____/____/ \\____/___//____/\n",
"\n",
"Welcome to GRASS GIS 7.2.1\n",
"GRASS GIS homepage: http://grass.osgeo.org\n",
"This version running through: Bash Shell (/bin/bash)\n",
"Help is available with the command: g.manual -i\n",
"See the licence terms with: g.version -c\n",
"See citation options with: g.version -x\n",
"Start the GUI with: g.gui wxpython\n",
"When ready to quit enter: exit\n",
"\n",
"Cleaning up temporary files...\n",
"Done.\n",
"\n",
"Goodbye from GRASS GIS\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"## Create a new location with geographic lat/lon grid\n",
"grass72 -c EPSG:4326 ~/Grass/wise\n",
"\n",
"## set the spatial extent and resulotion; save as default\n",
"g.region -sp s=-56 n=83 w=-180 e=180 res=0.5"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Data processing\n",
"\n",
"- Importing Arc/Info raster using GDAL (http://www.gdal.org/)\n",
"- conversion raster to vector polygon"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[H\u001b[2J +----------------------------------------------------------------------------+\n",
" | Map: wise_basemap Date: Tue Sep 5 09:35:16 2017 |\n",
" | Mapset: work Login of Creator: jsteinkamp |\n",
" | Location: wise |\n",
" | DataBase: /Users/jsteinkamp/Grass |\n",
" | Title: |\n",
" | Timestamp: none |\n",
" |----------------------------------------------------------------------------|\n",
" | |\n",
" | Type of Map: raster Number of Categories: 0 |\n",
" | Data Type: CELL |\n",
" | Rows: 360 |\n",
" | Columns: 720 |\n",
" | Total Cells: 259200 |\n",
" | Projection: Latitude-Longitude |\n",
" | N: 90N S: 90S Res: 0:30 |\n",
" | E: 180E W: 180W Res: 0:30 |\n",
" | Range of data: min = 1 max = 45948 |\n",
" | |\n",
" | Data Description: |\n",
" | generated by r.in.gdal |\n",
" | |\n",
" | Comments: |\n",
" | r.in.gdal --overwrite --quiet -o input=\"/Volumes/FAT/GIS/WISE_v3/Wis\\ |\n",
" | esnum/hdr.adf\" output=\"wise_basemap\" memory=300 offset=0 num_digits=0 |\n",
" | |\n",
" +----------------------------------------------------------------------------+\n",
"\n",
" +----------------------------------------------------------------------------+\n",
" | Name: wise |\n",
" | Mapset: work |\n",
" | Location: wise |\n",
" | Database: /Users/jsteinkamp/Grass |\n",
" | Title: |\n",
" | Map scale: 1:1 |\n",
" | Name of creator: jsteinkamp |\n",
" | Organization: |\n",
" | Source date: Tue Sep 5 09:35:16 2017 |\n",
" | Timestamp (first layer): none |\n",
" |----------------------------------------------------------------------------|\n",
" | Map format: native |\n",
" |----------------------------------------------------------------------------|\n",
" | Type of map: vector (level: 2) |\n",
" | |\n",
" | Number of points: 0 Number of centroids: 55690 |\n",
" | Number of lines: 0 Number of boundaries: 116792 |\n",
" | Number of areas: 55690 Number of islands: 132 |\n",
" | |\n",
" | Map is 3D: No |\n",
" | Number of dblinks: 1 |\n",
" | |\n",
" | Projection: Latitude-Longitude |\n",
" | |\n",
" | N: 83N S: 56S |\n",
" | E: 180E W: 180W |\n",
" | |\n",
" | Digitization threshold: 0 |\n",
" | Comment: |\n",
" | |\n",
" +----------------------------------------------------------------------------+\n",
"\n",
"\u001b[H\u001b[2J"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cleaning up temporary files...\n",
"Creating new GRASS GIS location/mapset...\n",
"Missing WIND file fixed\n",
"\n",
" __________ ___ __________ _______________\n",
" / ____/ __ \\/ | / ___/ ___/ / ____/ _/ ___/\n",
" / / __/ /_/ / /| | \\__ \\\\_ \\ / / __ / / \\__ \\\n",
" / /_/ / _, _/ ___ |___/ /__/ / / /_/ // / ___/ /\n",
" \\____/_/ |_/_/ |_/____/____/ \\____/___//____/\n",
"\n",
"Welcome to GRASS GIS 7.2.1\n",
"GRASS GIS homepage: http://grass.osgeo.org\n",
"This version running through: Bash Shell (/bin/bash)\n",
"Help is available with the command: g.manual -i\n",
"See the licence terms with: g.version -c\n",
"See citation options with: g.version -x\n",
"Start the GUI with: g.gui wxpython\n",
"When ready to quit enter: exit\n",
"\n",
"WARNING: Over-riding projection check\n",
"WARNING: File <map.png> already exists and will be overwritten\n",
"Output file: /Users/jsteinkamp/Grass/map.png\n",
"Cleaning up temporary files...\n",
"Done.\n",
"\n",
"Goodbye from GRASS GIS\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"grass72 -c ~/Grass/wise/work\n",
"\n",
"## set the path variable to the data files\n",
"INPUTDIR=\"/Volumes/FAT/GIS\"\n",
"\n",
"## import the Raster basemap from an ARC/Info binary grid \n",
"r.in.gdal -o --overwrite --quiet input=${INPUTDIR}/WISE_v3/Wisesnum/hdr.adf output=wise_basemap\n",
"r.info wise_basemap\n",
"\n",
"## convert the raster data to a polygon (vector)\n",
"r.to.vect --overwrite --quiet input=wise_basemap output=wise type=area\n",
"v.info wise\n",
"\n",
"## export to PNG\n",
"d.mon --overwrite start=png\n",
"d.rast --quiet wise_basemap\n",
"d.mon stop=png"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Importing the database\n",
"\n",
"Using commands from mdbtools to\n",
"- initialize the sqlite database table\n",
"- export to csv and import to sqlite"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[H\u001b[2JSNUM|AREA1|DEPT_1|AREA2|DEPT_2|AREA3|DEPT_3|AREA4|DEPT_4|AREA5|DEPT_5|AREA6|DEPT_6|AREA7|DEPT_7|AREA8|DEPT_8|AREA9|DEPT_9|AREA10|DEPT_10\n",
"1|15.3000002|145.0|14.8999996|120.0|14.1000004|10.0|13.8000002|120.0|13.8000002|117.0|7.19999981|116.0|6.0999999|150.0|6.0999999|140.0|5.5999999|120.0|3.0999999|130.0\n",
"\u001b[H\u001b[2J"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cleaning up temporary files...\n",
"Starting GRASS GIS...\n",
"\n",
" __________ ___ __________ _______________\n",
" / ____/ __ \\/ | / ___/ ___/ / ____/ _/ ___/\n",
" / / __/ /_/ / /| | \\__ \\\\_ \\ / / __ / / \\__ \\\n",
" / /_/ / _, _/ ___ |___/ /__/ / / /_/ // / ___/ /\n",
" \\____/_/ |_/_/ |_/____/____/ \\____/___//____/\n",
"\n",
"Welcome to GRASS GIS 7.2.1\n",
"GRASS GIS homepage: http://grass.osgeo.org\n",
"This version running through: Bash Shell (/bin/bash)\n",
"Help is available with the command: g.manual -i\n",
"See the licence terms with: g.version -c\n",
"See citation options with: g.version -x\n",
"Start the GUI with: g.gui wxpython\n",
"When ready to quit enter: exit\n",
"\n",
"No MSysRelationships\n",
"Cleaning up temporary files...\n",
"Done.\n",
"\n",
"Goodbye from GRASS GIS\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"grass72 ~/Grass/wise/work\n",
"\n",
"## set some variables to save typing\n",
"INPUTDIR=\"/Volumes/FAT/GIS\"\n",
"GrassDB=$(g.gisenv GISDBASE)/$(g.gisenv LOCATION_NAME)/$(g.gisenv MAPSET)/sqlite/sqlite.db\n",
"\n",
"## table layout export/import\n",
"mdb-schema ${INPUTDIR}/WISE_v3/WISE30by30min.mdb -T yDEPT sqlite | sqlite3 $GrassDB\n",
"\n",
"## table data export to CSV\n",
"mdb-export -d \"|\" ${INPUTDIR}/WISE_v3/WISE30by30min.mdb yDEPT > ${TMPDIR}/yDEPT.csv\n",
"\n",
"## table import\n",
"echo \".import ${TMPDIR}/yDEPT.csv yDEPT\" | sqlite3 $GrassDB\n",
"rm ${TMPDIR}/yDEPT.csv\n",
"echo -e \".header on\\n select * from yDEPT where SNUM==1;\" | sqlite3 $GrassDB"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## SQL joins and raster calculations\n",
"\n",
"The table layout, e.g. column names must be known in advance.\n",
"- Adding new columns to the georeferenced table\n",
"- filling it with the "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[H\u001b[2J\u001b[H\u001b[2J"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cleaning up temporary files...\n",
"Starting GRASS GIS...\n",
"\n",
" __________ ___ __________ _______________\n",
" / ____/ __ \\/ | / ___/ ___/ / ____/ _/ ___/\n",
" / / __/ /_/ / /| | \\__ \\\\_ \\ / / __ / / \\__ \\\n",
" / /_/ / _, _/ ___ |___/ /__/ / / /_/ // / ___/ /\n",
" \\____/_/ |_/_/ |_/____/____/ \\____/___//____/\n",
"\n",
"Welcome to GRASS GIS 7.2.1\n",
"GRASS GIS homepage: http://grass.osgeo.org\n",
"This version running through: Bash Shell (/bin/bash)\n",
"Help is available with the command: g.manual -i\n",
"See the licence terms with: g.version -c\n",
"See citation options with: g.version -x\n",
"Start the GUI with: g.gui wxpython\n",
"When ready to quit enter: exit\n",
"\n",
"Removing vector <wise>\n",
"Cleaning up temporary files...\n",
"Done.\n",
"\n",
"Goodbye from GRASS GIS\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"grass72 ~/Grass/wise/work\n",
"\n",
"## Set the path variable\n",
"GrassDB=$(g.gisenv GISDBASE)/$(g.gisenv LOCATION_NAME)/$(g.gisenv MAPSET)/sqlite/sqlite.db\n",
"\n",
"## for the sake of speed take only the three first ones, instead of all 10\n",
"## Ideally filter for real soils (values -95/-98/-99; rock/glacier/water)\n",
"for n in 1 2 3; do\n",
" sqlite3 ${GrassDB} \"ALTER TABLE wise ADD d${n} REAL\"\n",
" sqlite3 ${GrassDB} \"ALTER TABLE wise ADD a${n} REAL\"\n",
" sqlite3 ${GrassDB} \"UPDATE wise SET d${n}=(SELECT DEPT_${n} FROM yDEPT WHERE yDEPT.SNUM=wise.value)\"\n",
" sqlite3 ${GrassDB} \"UPDATE wise SET a${n}=(SELECT AREA${n} FROM yDEPT WHERE yDEPT.SNUM=wise.value)\"\n",
" v.to.rast --overwrite --quiet input=wise output=d${n} use=attr type=area attribute_column=d${n}\n",
" v.to.rast --overwrite --quiet input=wise output=a${n} use=attr type=area attribute_column=a${n}\n",
"\n",
" r.mapcalc --overwrite --quiet \"d${n} = if(isnull(d${n}), 0, d${n})\" \n",
" r.mapcalc --overwrite --quiet \"a${n} = if(isnull(a${n}), 0, a${n})\"\n",
" r.mapcalc --overwrite --quiet \"d${n} = if(d${n} < 0, 0, d${n})\"\n",
" r.mapcalc --overwrite --quiet \"a${n} = if(a${n} < 0, 0, a${n})\"\n",
"done\n",
"\n",
"## weighted mean; min/max\n",
"r.mapcalc --overwrite --quiet \"wise_coverage = a1 + a2 + a3\"\n",
"r.mapcalc --overwrite --quiet \"wise_depth = (d1*a1 + d2*a2 + d3*a3) / (a1 + a2 + a3)\"\n",
"r.mapcalc --overwrite --quiet \"wise_coverage = if(wise_coverage == 0, null(), wise_coverage)\"\n",
"r.mapcalc --overwrite --quiet \"wise_depth = if(wise_depth == 0, null(), wise_depth)\"\n",
"\n",
"## deleting old files\n",
"g.remove -f type=vect name=wise,yDEPT\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Exporting"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[H\u001b[2J\u001b[H\u001b[2J"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cleaning up temporary files...\n",
"Starting GRASS GIS...\n",
"\n",
" __________ ___ __________ _______________\n",
" / ____/ __ \\/ | / ___/ ___/ / ____/ _/ ___/\n",
" / / __/ /_/ / /| | \\__ \\\\_ \\ / / __ / / \\__ \\\n",
" / /_/ / _, _/ ___ |___/ /__/ / / /_/ // / ___/ /\n",
" \\____/_/ |_/_/ |_/____/____/ \\____/___//____/\n",
"\n",
"Welcome to GRASS GIS 7.2.1\n",
"GRASS GIS homepage: http://grass.osgeo.org\n",
"This version running through: Bash Shell (/bin/bash)\n",
"Help is available with the command: g.manual -i\n",
"See the licence terms with: g.version -c\n",
"See citation options with: g.version -x\n",
"Start the GUI with: g.gui wxpython\n",
"When ready to quit enter: exit\n",
"\n",
"Adding raster map <wise_coverage@work> to group\n",
"Adding raster map <wise_depth@work> to group\n",
"ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.\n",
"ERROR 6: SetColorTable() can only be called on band 1.\n",
"Warning 1: Metadata exceeding 32000 bytes cannot be written into GeoTIFF. Transferred to PAM instead.\n",
"Removing imagery group <soildepth>\n",
"Cleaning up temporary files...\n",
"Done.\n",
"\n",
"Goodbye from GRASS GIS\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"grass72 ~/Grass/wise/work\n",
"\n",
"OUTPUTDIR=$(g.gisenv GISDBASE)\n",
"\n",
"i.group group=soildepth input=wise_coverage,wise_depth\n",
"r.out.gdal --overwrite --quiet input=soildepth output=${OUTPUTDIR}/wise_soildepth.tif\n",
"g.remove -f type=group name=soildepth"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Display the source blob
Display the rendered blob
Raw
R-GIS.ipynb
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment