[GRASS-user] Simultaneous r.horizon processes

Daniel_Lee · April 1, 2012, 12:10pm

Hi list,

I’ve got a large area I’m running r.horizon on. Since the process takes quite a while but it only runs on one processor core and I’ve got plenty of RAM, I’ve decided to script it in Python and use run_command and start_command in order to create multiple horizon maps in parallel. Another reason I’m doing it that way is that I only need certain horizons, so I don’t want GRASS to do all 360°.

The problem is, when I run the script, all of the maps that I want made seem to be made and I can see in top that there are several instances of r.horizon. But only the instance to be finished last is actually saved. Here’s the output:

Calculating map 1 of 1 (angle 0.000000, raster map <horizon0_0>)
100%
Calculating map 1 of 1 (angle 12.500000, raster map <horizon12.5_0>)
Calculating map 1 of 1 (angle 7.500000, raster map <horizon7.5_0>)
Calculating map 1 of 1 (angle 42.500000, raster map <horizon42.5_0>)
Calculating map 1 of 1 (angle 150.000000, raster map <horizon150.0_0>)
100%
Calculating map 1 of 1 (angle 162.500000, raster map <horizon162.5_0>)
100%
Calculating map 1 of 1 (angle 175.000000, raster map <horizon175.0_0>)
Calculating map 1 of 1 (angle 167.500000, raster map <horizon167.5_0>)
Calculating map 1 of 1 (angle 197.500000, raster map <horizon197.5_0>)
Calculating map 1 of 1 (angle 205.000000, raster map <horizon205.0_0>)
Calculating map 1 of 1 (angle 225.000000, raster map <horizon225.0_0>)
100%
Calculating map 1 of 1 (angle 257.500000, raster map <horizon257.5_0>)
Calculating map 1 of 1 (angle 262.500000, raster map <horizon262.5_0>)
100%
Calculating map 1 of 1 (angle 275.000000, raster map <horizon275.0_0>)
100%
Calculating map 1 of 1 (angle 295.000000, raster map <horizon295.0_0>)
Calculating map 1 of 1 (angle 300.000000, raster map <horizon300.0_0>)
100%
Calculating map 1 of 1 (angle 322.500000, raster map <horizon322.5_0>)
Calculating map 1 of 1 (angle 325.000000, raster map <horizon325.0_0>)
100%
Calculating map 1 of 1 (angle 337.500000, raster map <horizon337.5_0>)
Calculating map 1 of 1 (angle 335.000000, raster map <horizon335.0_0>)
Calculating map 1 of 1 (angle 345.000000, raster map <horizon345.0_0>)
Calculating map 1 of 1 (angle 350.000000, raster map <horizon350.0_0>)
100%
Calculating map 1 of 1 (angle 352.500000, raster map <horizon352.5_0>)
100%
Calculating map 1 of 1 (angle 355.000000, raster map <horizon355.0_0>)
100%
Calculating map 1 of 1 (angle 357.500000, raster map <horizon357.5_0>)
100%
Calculating map 1 of 1 (angle 257.500000, raster map <horizon257.5_0>)
Calculating map 1 of 1 (angle 262.500000, raster map <horizon262.5_0>)
100%
Calculating map 1 of 1 (angle 275.000000, raster map <horizon275.0_0>)
100%
Calculating map 1 of 1 (angle 295.000000, raster map <horizon295.0_0>)
Calculating map 1 of 1 (angle 300.000000, raster map <horizon300.0_0>)
100%
Calculating map 1 of 1 (angle 322.500000, raster map <horizon322.5_0>)
Calculating map 1 of 1 (angle 325.000000, raster map <horizon325.0_0>)
100%
Calculating map 1 of 1 (angle 337.500000, raster map <horizon337.5_0>)
Calculating map 1 of 1 (angle 335.000000, raster map <horizon335.0_0>)
Calculating map 1 of 1 (angle 345.000000, raster map <horizon345.0_0>)
Calculating map 1 of 1 (angle 350.000000, raster map <horizon350.0_0>)
100%
Calculating map 1 of 1 (angle 352.500000, raster map <horizon352.5_0>)
100%
Calculating map 1 of 1 (angle 355.000000, raster map <horizon355.0_0>)
100%
Calculating map 1 of 1 (angle 357.500000, raster map <horizon357.5_0>)
100%

But when I do g.list rast at the end of it all I only have these maps:
raster files available in mapset :
dom_aspect horizon162.5_0 horizon325.0_0
dom_complete horizon225.0_0 horizon350.0_0
dom_slope horizon262.5_0 horizon352.5_0
horizon0_0 horizon275.0_0 horizon355.0_0
horizon150.0_0 horizon300.0_0 horizon357.5_0

As you can see, the only horizon maps saved are the last ones made. Does anybody know why this is the case and, if so, what would be a good workaround? Finding the horizons serially on only one processor means that my computer does a lot of sitting around twiddling thumbs. Thanks a bunch!

Daniel

–

B.Sc. Daniel Lee
Geschäftsführung für Forschung und Entwicklung
ISIS - International Solar Information Solutions GbR
Vertreten durch: Daniel Lee, Nepomuk Reinhard und Nils Räder

Deutschhausstr. 10
35037 Marburg
Festnetz: +49 6421 379 6256
Mobil: +49 176 6127 7269
E-Mail: Lee@isi-solutions.org
Web: http://www.isi-solutions.org

ISIS wird gefördert durch die Bundesrepublik Deutschland, Zuwendungsgeber: Bundesministerium für Wirtschaft und Technologie aufgrund eines Beschlusses des Deutschen Bundestages, sowie durch die Europäische Union, Zuwendungsgeber: Europäischer Sozialfonds.
Zusätzliche Unterstützung erhält ISIS von dem Entrepreneurship Cluster Mittelhessen, der Universität Marburg, dem Laboratory for Climatology and Remote Sensing und dem GIS-Lab Marburg.

Collin_Bode · April 17, 2012, 7:21pm

Lee,

I tried the same thing, trying to integrate a full workflow in Python for r.sun (see attached code). I am using a 2x2 meter lidar elevation map with about 120 million cells, so automation is pretty necessary. The whole thing takes 3 weeks to run using 8 CPU cores and 12GB ram.

However, r.horizon does not work properly on Grass 6.4.2. I can only get it to run on one core in sequence. I cannot get it to accept a single angle and run properly. The output is always for angle 0, both in name and output data. Bug? I was planning to post this to the dev list if no one responds here.

It would be worthwhile to find out if your version is simply overwriting the previous, similarly named map. If it is, and the maps are actually different (r.info), then you can change your code to run r.horizon for a single angle, then rename the map. Something like this:

for azimuth in range(0,361,step_size):
grass.run_command(“r.horizon”, elevin=dem,
direction=angle, maxdistance=maxdistance,
horizon=prefhor, dist=dist, overwrite = 1)
az = str(azimuth).zfill(3)
tempmap = prefhor+‘0’
horizonmap = prefhor+'’+az
cmd = tempmap+“,”+horizonmap
grass.run_command(“g.rename”,rast=cmd)

r.sun does work in multiple cores.

Collin Bode
Project Manager, Desktop Watershed Integrated Program
National Center for Earth-surface Dynamics
University of California, Berkeley

http://www.nced.umn.edu

http://ib.berkeley.edu/labs/power
http://angelo.berkeley.edu

(attachments)

2_rsun_mp.py (15.7 KB)

Collin_Bode · April 18, 2012, 4:54am

Daniel,

Your post prompted me to go back and retest the broken code I had for r.horizon. I have fixed it and have attached a script that I think is somewhat general, unlike the previous one.

There are several more issues I should mention:

Each process MUST run in its own mapset. Grass mapsets are not thread or process safe. The multiple instances will overwrite each other and create jibberish.
Use python scripts outside of grass so you can define the mapset when you create a separate process.
Manually spawn only enough processes as you have CPU cores available. After that, force your code to wait until that batch is completed. Once that set is done, spawn the next set. I did not do this in the code I posted earlier.

Multiprocessing with grass requires a lot of careful hand-holding to make it work.

hopefully the following script is useful.

Collin Bode
Project Manager, Desktop Watershed Integrated Program
National Center for Earth-surface Dynamics
University of California, Berkeley
Cell: 415-305-5346 Lab: 510-643-9294
http://www.nced.umn.edu

http://angelo.berkeley.edu
http://ncalm.berkeley.edu

(attachments)

horizon_mp.py (7.22 KB)

neteler · April 18, 2012, 10:21am

On Wed, Apr 18, 2012 at 6:54 AM, Collin Bode <collin@berkeley.edu> wrote:
...

There are several more issues I should mention:
1. Each process MUST run in its own mapset. Grass mapsets are not thread
or process safe. The multiple instances will overwrite each other and
create jibberish.

Right - same thing on a cluster. The solution is easy, see
"# generate machine (blade) unique TMP string"
in
http://grass.osgeo.org/wiki/Parallel_GRASS_jobs#The_Grid_Engine_script

which I have just updated (simplified new approach).

Markus

Daniel_Lee · April 18, 2012, 6:44pm

Hi Collin, Markus,

Thanks a lot for the script and tips, I really appreciate it. That’s definitely a much bigger chunk to chew on than I’d planned, but it’s good to know that. It seems like most of the stuff would be fairly simple to implement - the program I’m working with is built modularly. I’ve got a couple of questions that would give me some clarity here on the way though, if you don’t mind:

I notice in your script that you use the multiprocessing library? Are there any advantages in doing so? It looks to me like you’d be able to accomplish the same things without using the library so I’d be interested to know why you chose it.

What do you mean, r.sun can do multithreading? I’ve heard that r.sun uses multithreading in GRASS 7, but is that implemented in GRASS 6? Or are you talking about “poor man’s multithreading,” like on the GRASS wiki?

I haven’t noticed any problems with running r.horizon on GRASS 6.4.2, the angles looked correct in test data I’ve used thus far. The maps are all assigned different names later on and the results look fine, but a lot of them were missing. I was puzzled for a while because test data worked and my actual data didn’t, but I did a little debugging today and think it’s due to the large size of my study area (~36 km², 0.5m resolution). If I spatially partition the area and then stitch everything back together, I hope it works - tests on smaller regions have worked correctly thus far but I’ll need to wait a while to see the real results.

To make a long story short, I’m going to try using the tricks in Markus’ link to get the processes to coexist peacefully and then partition the area spatially so that it doesn’t kill my memory. Hopefully I’ll get better results then. In the meantime you’ve been a great help

Best,
Daniel

–

B.Sc. Daniel Lee
Geschäftsführung für Forschung und Entwicklung
ISIS - International Solar Information Solutions GbR
Vertreten durch: Daniel Lee, Nepomuk Reinhard und Nils Räder

Softwarecenter 3
35037 Marburg
Festnetz: +49 6421 379 6256
Mobil: +49 176 6127 7269
E-Mail: Lee@isi-solutions.org
Web: http://www.isi-solutions.org

Am 17. April 2012 21:21 schrieb Collin Bode <collin@berkeley.edu>:

Lee,

I tried the same thing, trying to integrate a full workflow in Python for r.sun (see attached code). I am using a 2x2 meter lidar elevation map with about 120 million cells, so automation is pretty necessary. The whole thing takes 3 weeks to run using 8 CPU cores and 12GB ram.

However, r.horizon does not work properly on Grass 6.4.2. I can only get it to run on one core in sequence. I cannot get it to accept a single angle and run properly. The output is always for angle 0, both in name and output data. Bug? I was planning to post this to the dev list if no one responds here.

It would be worthwhile to find out if your version is simply overwriting the previous, similarly named map. If it is, and the maps are actually different (r.info), then you can change your code to run r.horizon for a single angle, then rename the map. Something like this:

for azimuth in range(0,361,step_size):
grass.run_command(“r.horizon”, elevin=dem,
direction=angle, maxdistance=maxdistance,
horizon=prefhor, dist=dist, overwrite = 1)
az = str(azimuth).zfill(3)
tempmap = prefhor+‘0’
horizonmap = prefhor+'’+az
cmd = tempmap+“,”+horizonmap
grass.run_command(“g.rename”,rast=cmd)

r.sun does work in multiple cores.

Collin Bode
Project Manager, Desktop Watershed Integrated Program
National Center for Earth-surface Dynamics
University of California, Berkeley

http://www.nced.umn.edu

http://ib.berkeley.edu/labs/power
http://angelo.berkeley.edu

On Apr 1, 2012, at 5:10 AM, Daniel Lee wrote:

Hi list,

I’ve got a large area I’m running r.horizon on. Since the process takes quite a while but it only runs on one processor core and I’ve got plenty of RAM, I’ve decided to script it in Python and use run_command and start_command in order to create multiple horizon maps in parallel. Another reason I’m doing it that way is that I only need certain horizons, so I don’t want GRASS to do all 360°.

The problem is, when I run the script, all of the maps that I want made seem to be made and I can see in top that there are several instances of r.horizon. But only the instance to be finished last is actually saved. Here’s the output:

Calculating map 1 of 1 (angle 0.000000, raster map <horizon0_0>)
100%
Calculating map 1 of 1 (angle 12.500000, raster map <horizon12.5_0>)
Calculating map 1 of 1 (angle 7.500000, raster map <horizon7.5_0>)
Calculating map 1 of 1 (angle 42.500000, raster map <horizon42.5_0>)
Calculating map 1 of 1 (angle 150.000000, raster map <horizon150.0_0>)
100%
Calculating map 1 of 1 (angle 162.500000, raster map <horizon162.5_0>)
100%
Calculating map 1 of 1 (angle 175.000000, raster map <horizon175.0_0>)
Calculating map 1 of 1 (angle 167.500000, raster map <horizon167.5_0>)
Calculating map 1 of 1 (angle 197.500000, raster map <horizon197.5_0>)
Calculating map 1 of 1 (angle 205.000000, raster map <horizon205.0_0>)
Calculating map 1 of 1 (angle 225.000000, raster map <horizon225.0_0>)
100%
Calculating map 1 of 1 (angle 257.500000, raster map <horizon257.5_0>)
Calculating map 1 of 1 (angle 262.500000, raster map <horizon262.5_0>)
100%
Calculating map 1 of 1 (angle 275.000000, raster map <horizon275.0_0>)
100%
Calculating map 1 of 1 (angle 295.000000, raster map <horizon295.0_0>)
Calculating map 1 of 1 (angle 300.000000, raster map <horizon300.0_0>)
100%
Calculating map 1 of 1 (angle 322.500000, raster map <horizon322.5_0>)
Calculating map 1 of 1 (angle 325.000000, raster map <horizon325.0_0>)
100%
Calculating map 1 of 1 (angle 337.500000, raster map <horizon337.5_0>)
Calculating map 1 of 1 (angle 335.000000, raster map <horizon335.0_0>)
Calculating map 1 of 1 (angle 345.000000, raster map <horizon345.0_0>)
Calculating map 1 of 1 (angle 350.000000, raster map <horizon350.0_0>)
100%
Calculating map 1 of 1 (angle 352.500000, raster map <horizon352.5_0>)
100%
Calculating map 1 of 1 (angle 355.000000, raster map <horizon355.0_0>)
100%
Calculating map 1 of 1 (angle 357.500000, raster map <horizon357.5_0>)
100%
Calculating map 1 of 1 (angle 257.500000, raster map <horizon257.5_0>)
Calculating map 1 of 1 (angle 262.500000, raster map <horizon262.5_0>)
100%
Calculating map 1 of 1 (angle 275.000000, raster map <horizon275.0_0>)
100%
Calculating map 1 of 1 (angle 295.000000, raster map <horizon295.0_0>)
Calculating map 1 of 1 (angle 300.000000, raster map <horizon300.0_0>)
100%
Calculating map 1 of 1 (angle 322.500000, raster map <horizon322.5_0>)
Calculating map 1 of 1 (angle 325.000000, raster map <horizon325.0_0>)
100%
Calculating map 1 of 1 (angle 337.500000, raster map <horizon337.5_0>)
Calculating map 1 of 1 (angle 335.000000, raster map <horizon335.0_0>)
Calculating map 1 of 1 (angle 345.000000, raster map <horizon345.0_0>)
Calculating map 1 of 1 (angle 350.000000, raster map <horizon350.0_0>)
100%
Calculating map 1 of 1 (angle 352.500000, raster map <horizon352.5_0>)
100%
Calculating map 1 of 1 (angle 355.000000, raster map <horizon355.0_0>)
100%
Calculating map 1 of 1 (angle 357.500000, raster map <horizon357.5_0>)
100%

But when I do g.list rast at the end of it all I only have these maps:
raster files available in mapset :
dom_aspect horizon162.5_0 horizon325.0_0
dom_complete horizon225.0_0 horizon350.0_0
dom_slope horizon262.5_0 horizon352.5_0
horizon0_0 horizon275.0_0 horizon355.0_0
horizon150.0_0 horizon300.0_0 horizon357.5_0

As you can see, the only horizon maps saved are the last ones made. Does anybody know why this is the case and, if so, what would be a good workaround? Finding the horizons serially on only one processor means that my computer does a lot of sitting around twiddling thumbs. Thanks a bunch!

Daniel

–

B.Sc. Daniel Lee
Geschäftsführung für Forschung und Entwicklung
ISIS - International Solar Information Solutions GbR
Vertreten durch: Daniel Lee, Nepomuk Reinhard und Nils Räder

Deutschhausstr. 10
35037 Marburg
Festnetz: +49 6421 379 6256
Mobil: +49 176 6127 7269
E-Mail: Lee@isi-solutions.org
Web: http://www.isi-solutions.org

ISIS wird gefördert durch die Bundesrepublik Deutschland, Zuwendungsgeber: Bundesministerium für Wirtschaft und Technologie aufgrund eines Beschlusses des Deutschen Bundestages, sowie durch die Europäische Union, Zuwendungsgeber: Europäischer Sozialfonds.
Zusätzliche Unterstützung erhält ISIS von dem Entrepreneurship Cluster Mittelhessen, der Universität Marburg, dem Laboratory for Climatology and Remote Sensing und dem GIS-Lab Marburg.

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user