[GRASS-user] Best workflow

Paul_Meems2 · June 24, 2016, 7:13am

Hi list,

I need to do the following and I want to know what is the best workflow.

I have a PostGIS database with polygons.
I have a folder with height-maps in tiff-format.
For each polygon I need to get the height-map.

This smaller height map needs to send to GRASS and several commands need to be performed.
At the end GRASS will produce a vector file. This vector data needs to be saved back in PostGIS.

I already did all the steps manually but need to automate it because I will process all (7M+) polygons. I was planning to set up an Ubuntu server with the PostGIS database.
Because it is possible a polygon is on the edge of multiple tiff-files I thought of setting up GeoServer and feed it all tif-files and use WCS to get the needed height map.
Next do the GRASS commands. I want to use Python to perform all actions.

I’m rather new to GRASS. Is the above workflow the best approach? We will do all actions once for all polygons. But performance is an issue.

Hope to get any pointers or advice.

Thanks,

···

Paul Meems

The Netherlands

Blumentrath_Stefan1 · June 24, 2016, 8:17am

From: grass-user [mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Paul Meems

···

Depends on several things, e.g

do polygons overlap or not,
do “height maps” overlap or not?
What are the “several commands” you need to be performed

On the basis of the information you provided I would have a look at

gdalbuildvrt for merging your “height map” tiles
import PostGIS polygon table to GRASS
convert to raster using v.to.rast and assign ID to each raster representation of your polygons
calculate statistics using either

o r.stats

o r.univar

to create CSV files, which can be read back to PostGIS (COPY command) and then joined to your original map using ID field

I would not consider WCS / GeoServer, that gives you only overhead and extra work…

Cheers

Stefan

Paul_Meems2 · June 24, 2016, 9:34am

Thanks Stefan,

I didn’t thought of gdalbuildvrt. I’ll try that option.
The polygons do not overlap and I’m not sure about the height data. Perhaps they overlap a bit.

What I need to do in GRASS is:

call r.slope.aspect to calculate the slope and aspect. I had planned doing this for every polygon but when I’m going to use gdalbuildvrt I can do it for the whole dataset in one time.
Is that a good choice? The height data spans The Netherlands and has 0.5m resolution. Can GRASS handle such a big dataset?
Next command is r.mapcalc. This can be done on a subset or on the whole dataset at once.
Then I need to call r.to.vect. I think it is best to do this for a subset because I only need the part within the polygon. I think I can set the region to the polygon and then call r.to.vect to speed this up, right?
Then I need to clip with the border and export this vector dataset to PostGIS.
With the above workflow I don’t think I need to convert the polygons to raster first, right?

I’m reading https://grasswiki.osgeo.org/wiki/PostGIS and I’m using GRASS v7 so if I understand it correctly I can link a PostGIS table using v.external
Can I then loop through each feature and perform the above workflow and also update a field in PostGIS called something like ‘processed’ and set it to True?And am I right I can do this all using Python? I’m running Ubuntu server so no GUI available.

Thanks

···

2016-06-24 10:17 GMT+02:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

From: grass-user [mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Paul Meems

Depends on several things, e.g

do polygons overlap or not,

do “height maps” overlap or not?

What are the “several commands” you need to be performed

On the basis of the information you provided I would have a look at

gdalbuildvrt for merging your “height map” tiles

import PostGIS polygon table to GRASS

convert to raster using v.to.rast and assign ID to each raster representation of your polygons

calculate statistics using either

o r.stats

o r.univar

to create CSV files, which can be read back to PostGIS (COPY command) and then joined to your original map using ID field

I would not consider WCS / GeoServer, that gives you only overhead and extra work…

Cheers

Stefan

Paul

Blumentrath_Stefan1 · June 24, 2016, 9:55am

Hi Paul,

Sure, GRASS and Python is a powerful combination. See: https://grasswiki.osgeo.org/wiki/GRASS_and_Python and https://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library

But there are others on this list that are much more competent than I am on this topic.

GRASS 7 is definitely the right choice for heavy datasets.

However, is it really necessary to use 0.5m resolution for you purpose?

If yes, have a look at this: https://grass.osgeo.org/grass70/manuals/libpython/pygrass.modules.grid.html

Or this: https://grasswiki.osgeo.org/wiki/Parallel_GRASS_jobs

Not sure about the performance for r.external and VRT mosaics… So, you might be better off running parallel jobs for each tile (e.g. each tile in it`s own mapset)… And that you finally collect the results for all tiles…

If you want to assign raster values to polygons you will have to rasterize them (which is v.rast.stats https://grass.osgeo.org/grass70/manuals/v.rast.stats.html does).

Keep in mind that v.external does not build a proper topology. That might be unsafe…

Cheers

Stefan

···

Thanks Stefan,

I didn’t thought of gdalbuildvrt. I’ll try that option.

The polygons do not overlap and I’m not sure about the height data. Perhaps they overlap a bit.

What I need to do in GRASS is:

call r.slope.aspect to calculate the slope and aspect. I had planned doing this for every polygon but when I’m going to use gdalbuildvrt I can do it for the whole dataset in one time.
Is that a good choice? The height data spans The Netherlands and has 0.5m resolution. Can GRASS handle such a big dataset?
Next command is r.mapcalc. This can be done on a subset or on the whole dataset at once.
Then I need to call r.to.vect. I think it is best to do this for a subset because I only need the part within the polygon. I think I can set the region to the polygon and then call r.to.vect to speed this up, right?
Then I need to clip with the border and export this vector dataset to PostGIS.

With the above workflow I don’t think I need to convert the polygons to raster first, right?

I’m reading https://grasswiki.osgeo.org/wiki/PostGIS and I’m using GRASS v7 so if I understand it correctly I can link a PostGIS table using v.external

Can I then loop through each feature and perform the above workflow and also update a field in PostGIS called something like ‘processed’ and set it to True?

And am I right I can do this all using Python? I’m running Ubuntu server so no GUI available.

Thanks

Paul

2016-06-24 10:17 GMT+02:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

From: grass-user [mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Paul Meems

Depends on several things, e.g

do polygons overlap or not,

do “height maps” overlap or not?

What are the “several commands” you need to be performed

On the basis of the information you provided I would have a look at

gdalbuildvrt for merging your “height map” tiles

import PostGIS polygon table to GRASS

convert to raster using v.to.rast and assign ID to each raster representation of your polygons

calculate statistics using either

o r.stats

o r.univar

to create CSV files, which can be read back to PostGIS (COPY command) and then joined to your original map using ID field

I would not consider WCS / GeoServer, that gives you only overhead and extra work…

Cheers

Stefan

Blumentrath_Stefan1 · June 24, 2016, 10:44am

Forgot to mention: pay attention to border effects in r.slope.aspect (see: https://grass.osgeo.org/grass70/manuals/r.slope.aspect.html). So if your tiles do not overlap, you might have to merge them…

neteler · June 25, 2016, 11:34pm

On Fri, Jun 24, 2016 at 12:44 PM, Blumentrath, Stefan
<Stefan.Blumentrath@nina.no> wrote:

Forgot to mention: pay attention to border effects in r.slope.aspect (see:
https://grass.osgeo.org/grass70/manuals/r.slope.aspect.html). So if your
tiles do not overlap, you might have to merge them...

If you can define your tiles yourself, then you may use the "overlap"
parameter to create tiles which overlap by 1 row/column:

https://grass.osgeo.org/grass72/manuals/r.tile.html

Markus

Paul_Meems2 · June 28, 2016, 6:32am

Thanks Markus for your suggestion.

I first tried Stefan’s suggestion of using a VRT-file.
I have 1300 tiff-files with a combined file size of 400GB. Creating the VRT-file was no problem, but when I tried to import the VRT-file in GRASS7 using r.in.gdal GRASS crashes and dies.
Loading the same VRT-file in QGis did take some time to load but it worked.

So I think I need to set-up a GeoServer after all, just to get a small portion of my height data to process in GRASS without worrying about polygons that are partially on multiple tiff-files.

Of course I’m open to any other suggestion. As said I’m new to GRASS.

···

2016-06-26 1:34 GMT+02:00 Markus Neteler <neteler@osgeo.org>:

On Fri, Jun 24, 2016 at 12:44 PM, Blumentrath, Stefan
<Stefan.Blumentrath@nina.no> wrote:

Forgot to mention: pay attention to border effects in r.slope.aspect (see:
https://grass.osgeo.org/grass70/manuals/r.slope.aspect.html). So if your
tiles do not overlap, you might have to merge them…

If you can define your tiles yourself, then you may use the “overlap”
parameter to create tiles which overlap by 1 row/column:

https://grass.osgeo.org/grass72/manuals/r.tile.html

Markus

Paul

Paul Meems
Release manager, configuration manager
and forum moderator of MapWindow GIS.
www.mapwindow.org

Owner of MapWindow.nl - Support for
Dutch speaking users.
www.mapwindow.nl

Download the latest MapWinGIS mapping engine.

Download the latest MapWindow 5 open source desktop application.

neteler · June 28, 2016, 9:04pm

Hi Paul,

On Tue, Jun 28, 2016 at 8:32 AM, Paul Meems <bontepaarden@gmail.com> wrote:

Thanks Markus for your suggestion.

I first tried Stefan's suggestion of using a VRT-file.
I have 1300 tiff-files with a combined file size of 400GB. Creating the VRT-file was no problem, but when I tried to import the VRT-file in GRASS7 using r.in.gdal GRASS crashes and dies.

Huh - should not happen .. But which operating system are you using,
which GRASS version etc?
You may post the result of

g.version -bge

Loading the same VRT-file in QGis did take some time to load but it worked.

Please post also the output of

gdalinfo yourfile.vrt

So I think I need to set-up a GeoServer after all, just to get a small portion of my height data to process in GRASS without worrying about polygons that are partially on multiple tiff-files.

Of course I'm open to any other suggestion. As said I'm new to GRASS.

Let's try harder, data in the your file size range should not be a problem.

Markus