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@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