I am using GRASS r.skyview command (v.7, Windows) to determine "urban canyon effect" on pollution due to the presence of buildings. The command nicely works, unfortunately I need to process a very large raster (representing buildings in an entire municipality) with a high resolution (10260x6762 pixels) and the command, which process every single pixel of the raster, run very slowly (about 24h to run on the entire area).

Since in the end I need the value of the skyview factor only for some specific points in space, I am wondering if it is possible to automatically do de following:

1) import in GRASS a vector layer representing point locations

2) select a single point

3) set GRASS region to a 100m x 100m square around selected point

4) launch r.skyview command

and loop through all points in vector layer.

  • Alternative approaches are discussed in the thread at gis.stackexchange.com/a/7655. – whuber Feb 20 '14 at 17:48
  • Did you try the grass region iteration in my answer already? I'm curious about if this could be an implementation for large raster with high resolutions. – Stefan Mar 4 '14 at 9:01
  • Thank you @StefanB. for your interest. I resolved with a Python script (for Windows) that iterate only the definition of the region through the command g.region . First I import point coordinates as a list throw NumPy . Then I constructed a cycle, iterating on each point, that (i) define the GRASS region through the command g.region, setting parameters n=(point.Ycoord + 100), s=(point.Xcoord - 100) etc., (ii) run the command r.skyview in the defined region only, (iii) save the output map as a new raster. – mic_cord Mar 5 '14 at 9:15
  • @mic_cord. It would be helpful to dilate your short description here to an answer, even if you are answering your question yourself. You spoke about a square region. Are there some "effects" or differences between the use of square or round regions for your implementation? – Stefan Mar 5 '14 at 9:46

I've found a solution for LINUX. The bash scripting here is very comfortable. So you have to transfer the ideas explained here to windows. For this purpose you can use cygwin.

  1. Import your vector layer with the point features:

    v.in.ogr -o dsn=/your_shapefile_directory layer=your_vector_layer output=points
  2. The same for your raster map. I've imported an elevation model (25m DEM) for testing. I'm not using GRASS GIS 7 and cannot use r.skyview, so I calculated slope and aspect cells.

    r.in.gdal input=/.../...dem.tif output=dem
  3. You can create a 100x100 square region with v.mkgrid. Before you have to move your points with half the value of your square region, because the point is aligned to the lower left corner of 100x100m square (rectangle). So you have to move your points with v.edit 50m to the west and 50m to the south. Besides you have to set the parameter ids, cats or where. I have an ID column that defines my points exactly and I defined the where clause where="(id > 1)".

    v.edit map=points type=point tool=move move=-50,-50 where="(id > 1)"
  4. Now, you can create the grids for every point in the layer that are needed for the regions

    for i in `v.out.ascii input=points format=point fs=,`;
    do k=$((k+1)); v.mkgrid map="pp$k" grid=1,1 position=coor coor=${i%,*} box=100,100; done

    The variable k is for the names of the vector layers which are containing the single grids. In my approach I have 526 points and there are 526 resulting grids. The grids are labeled pp1 to pp526. In the v.mkgrid function there you can define coordinates on which the grid will be snapped. To get the coordinates from the points I've used v.out.ascii. A sample output of this is 4462191,5453625,499(these are Gauß-Krueger coordinates). Defining coor=${i%,*} in v.mkgrid sets the coordinates for the grids (grid=1,1 > number of rows and columns in grid), avoiding the ,499. In this picture you see the moved point (blue), the original point (red) and the processed grid (black box). enter image description here

  5. Now my answer goes further with the calculation of a slope for every little region defined by the grid. You have to use here your r.skyview.

    for i in `g.mlist type=vect pattern="pp*"`; 
    do g.region vect=$i; r.slope.aspect elevation=dem slope="s$i" aspect="a$i"; done

    With g.mlist I call all grids for setting up the right region for the raster processing with r.slope.aspect. I have 526 little raster parcels which I have to merge (one slope parcel called e.g. spp2)

  6. Combining the raster parcels. It is important to set the region for every raster (slope) parcel.

    for i in `g.mlist type=rast pattern="spp*"`; do names="$i,$names_rast"; done
    g.region rast=${names_rast%,}
    r.patch --o input=${names_rast%,} output=merged_slope

The result is one raster file with the raster parcels (here a slope raster). The raster patching for 526 raster parcels takes around 20min. When you don't set the region for the raster patching, there appears a warning: WARNING: Not enough room in history file for command line (truncated).. You can combine the raster parcels with r.series --o input=${names_rast%,} output=merged_slope method=sum, as well.

enter image description here

Because of the 25m resolution,the 100x100m parcels are not filled out with raster data.

  • I've edited my answer. I forgot to set the region in step 6 and had to use an extra script, because of the warning. Now it is more simple. – Stefan Feb 21 '14 at 13:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.