[GRASS-user] Rasterize polygons, multiple polygons per grid cell

Micha Silver micha at arava.co.il
Sat Sep 29 08:09:32 PDT 2012


On 09/28/2012 09:49 PM, Lars Dalby wrote:
> Hi
>
> I have been desperately trying to figure out how to solve the task described
> below in GRASS (via R using spgrass6). I have searched the list, but did not
> find exactly what I was looking for.
>
> I have a large polygon file, app. 250000 polygons. It is the GLWD-2 data set
> available  here
> <http://worldwildlife.org/publications/global-lakes-and-wetlands-database-small-lake-polygons-level-2>
> The map is in lat long. My goal is to end up with a raster layer in Behrmann
> projection with app 96486.28m by 96486.28m grid cells showing the area per
> grid cell covered by polygons in the GLWD-2 layer.

If I understand your task, you want to
* import the GLWD lakes polygon shapefile into GRASS
* reproject to a projected CRS
* create a GRID with cell size about 92 km. X 92 km, and
* sum up the area of lakes within each grid cell.
Here's a GRASS-only procedure to do that:

First I downloaded the polygon shapefile you linked to. THen I started 
GRASS in a WGS84 based LOCATION, and used:
# Use the '-o' option since the downloaded shapefile does not have a *.prj
v.in.ogr -o dsn=glwd_2.shp out=glwd

Next I restarted GRASS and created a new location based on the EPSG code 
3925. This is NOT World Behrmann, but it's a cylindrical equal area 
projection with meters as units, like Behrmann. I couldn't find any 
proj4 reference to Behrmann :-(

At this step, I made sure to set the default database connection to 
sqlite. (This is important later to aggregate the water body areas in 
each grid cell.):
eval `g.gisenv`
db.connect driver=sqlite database=$GISDBASE/$LOCATION_NAME/$MAPSET/sqlite.db

Now:
v.proj in=glwd loc=<the WGS84 LOCATION> map=<the previous MAPSET>

Next step, I created a grid using v.mkgrid. I limited the region to just 
Africa so things would move along a little faster. So something like:
# Create a vector grid with size 92,000 meters by 92,000 meters
v.mkgrid --o grid92 position=coor coor=-2000000,-5000000 box=92000,92000 
grid=140,85
# I got the 'grid=' parameters by dividing the length and width by 92,000

Now I use v.overlay to find the intersection of these two vectors. This 
takes a while since all the lakes need to be clipped at the boundaries 
of the grid cells.
v.overlay ain=grid92 bin=glwd out=glwd_grid operator=and

After this finishes, I add a new column for the areas of each water body 
(now clipped to each grid cell):

v.db.addcol glwd_grid col="clip_area double precision"
# calculate the area of each clipped water surface
v.to.db glwd_grid opt=area col=clip_area

The new glwd_grid vector now has in its attrib table:
v.info -c glwd_grid
  Displaying column types/names for database connection of layer 1:
INTEGER|cat
INTEGER|a_cat
INTEGER|a_row
INTEGER|a_col
INTEGER|b_cat
INTEGER|b_GLWD_ID
CHARACTER|b_TYPE
CHARACTER|b_POLY_SRC
DOUBLE PRECISION|b_AREA_SKM
DOUBLE PRECISION|b_PERIM_KM
DOUBLE PRECISION|b_LONG_DEG
DOUBLE PRECISION|b_LAT_DEG
DOUBLE PRECISION|clip_area

What we need to do is SUM all clip_area values for each a_cat (grid cell 
category). I did this right in sqlite as follows:
sqlite3 ~/geodata/grass/Behrmann/PERMANENT/sqlite.db
SQLite version 3.6.20
Enter ".help" for instructions
Enter SQL statements terminated with a ";"
sqlite> SELECT a_row AS Grid_Row, a_col AS Grid_Col, SUM(clip_area) AS 
Sum_Areas FROM glwd_grid GROUP BY a_cat ORDER BY a_row, a_col;
Grid_Row    Grid_Col    Sum_Areas
----------  ----------  ---------------
9           43          11734990.811648
10          42          7600679.947002
10          43          20154750.823557
10          44          10095944.174832
10          46          18364623.744757
10          49          3054521.141835
11          41          250083.456892
11          42          9839134.964962
11          43          23519926.337785
11          44          11323614.28017
11          45          8426141.099425
11          46          12079532.459832
11          47          11135300.910967
11          49          571804.968053
11          50          1319664.199523
12          41          18036461.225015
12          42          15664405.025981
12          43          2471761.404269
12          45          2339714.549891
12          47          2890312.572311
12          48          5728571.755616
12          49          7189360.100137
12          50          3492556.03765
12          51          7318834.4634


Here's the same query using db.select:
echo 'SELECT a_row AS GridRow, a_col AS GridCol, SUM(clip_area) As 
"Water Area in Grid" FROM glwd_grid GROUP BY a_cat ORDER BY a_row, a_col 
LIMIT 20;' | db.select
GridRow|GridCol|Water Area in Grid
9|43|11734990.811648
10|42|7600679.947002
10|43|20154750.823557
10|44|10095944.174832
10|46|18364623.744757
10|49|3054521.141835
11|41|250083.456892
11|42|9839134.964962
11|43|23519926.337785
11|44|11323614.28017
11|45|8426141.099425
11|46|12079532.459832
11|47|11135300.910967
11|49|571804.968053
11|50|1319664.199523
12|41|18036461.225015
12|42|15664405.025981
12|43|2471761.404269
12|45|2339714.549891
12|47|2890312.572311


HTH,
Micha

> So as I see it the task is to, for each grid cell, "clip" the underlying
> polygon file, calculate the area and return that value in the resulting
> raster.
>
> So far I have successfully imported the shape file and reprojected it to a
> location with the desired projection and resolution. So far so good.
>
> Then I imported a raster with the desired grid cell size and dimensions. So
> now I seem to have the layers I need in the same location.
>
> Then I some how want to do the actual clipping and calculation of the area
> of polygons in each cell. This is where I get stuck. Any suggestions to a
> newbie?
>
> Initially I thought I somehow could specify which function(e.g. sum, mean
> etc) should be used in v.to.rast in the case of multiple polygons falling in
> the same grid cell. But I could not find anything here.
>
> I'm running GRASS 7.0 on OSX btw.
>
> Any help and suggestions would be highly appreciated!
>
> Best regards,
> Lars
>
>
>
> --
> View this message in context: http://osgeo-org.1560.n6.nabble.com/Rasterize-polygons-multiple-polygons-per-grid-cell-tp5005264.html
> Sent from the Grass - Users mailing list archive at Nabble.com.
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
>
> This mail was received via Mail-SeCure System.
>
>


-- 
Micha Silver
GIS Consultant, Arava Development Co.
http://www.surfaces.co.il



More information about the grass-user mailing list