[GRASS-user] Best workflow

[Blumentrath, Stefan]

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


From: Paul Meems [mailto:bontepaarden at gmail.com]
Sent: 24. juni 2016 11:35
To: GRASS user list <grass-user at lists.osgeo.org>
Cc: Blumentrath, Stefan <Stefan.Blumentrath at nina.no>
Subject: Re: [GRASS-user] Best workflow

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 at nina.no<mailto:Stefan.Blumentrath at nina.no>>:
> From: grass-user [mailto:grass-user-bounces at lists.osgeo.org<mailto:grass-user-bounces at 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20160624/6abf8579/attachment-0001.html>