You need to understand the difference between:
- grass.run_command which performs a GRASS command
grass.run_command('g.remove', raster="raster1")
- grass.read_command which is used to receive a result (generally a string, same result as the command in GRASS GIS, of little interest to Python, you need to use the re (Regular Expression) module to process the result)
grass.read_command("g.region",flags="p")
'projection: 99 (Lambert Conformal Conic)\nzone: 0\ndatum: bel72\nellipsoid: international\nnorth: 112990.07863053\nsouth: 108258.30027318\nwest: 227552.40388837\neast: 235027.97855494\nnsres: 6.35138035\newres: 6.35138035\nrows: 745\ncols: 1177\ncells: 876865\n'
- grass.parse_command that parse the result for Python (generally a dictionary)
grass.parse_command("g.region",flags="p")
{'ellipsoid: international': None, 'zone: 0': None, 'nsres: 6.35138035': None, 'south: 108258.30027318': None, 'ewres: 6.35138035': None, 'west: 227552.40388837': None, 'rows: 745': None, 'east: 235027.97855494': None, 'projection: 99 (Lambert Conformal Conic)': None, 'cells: 876865': None, 'north: 112990.07863053': None, 'cols: 1177': None, 'datum: bel72': None}
So, in your case, it is not:
grass.run_command("g.list", _type="rast")
whose result is 0
but:
rasters = grass.read_command("g.list",_type="rast")
whose result is a string with all the layers in the variable rasters
type(raster)
<type 'str'>
But it is not interesting at all here, because you only need the names of the layers:
rasters = grass.parse_command("g.list", _type="rast")
whose result is a Python dictionary and for further, see Python Scripts For GRASS GIS: Using "parse_command()"
The result is a list with the layer names:
['raster1', 'raster2', 'raster3']
so the processing is:
for items in mymaps:
grass.parse_command('r.univar', map=items)
.
Since you put GRASS into the tags, here's a solution based on GRASS:
First, you need to know exactly what coordinate system the original data is in (as always with GRASS). I see that the *.prj file contains "TRANSVERSE MERCATOR" but it's not one of the standard UTM zones. Since you have not mentioned the EPSG code of this data, here's the proj4 string for creating a new matching GRASS Location:
+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +units=m +no_defs
Now start GRASS, create a LOCATION (based on the above) and a MAPSET.
And now import the shapefile into GRASS by running:
v.in.ogr input="sp_data.shp" output=sp_data
Add a column for the area of each feature, and calculate its area:
v.db.addcolumn map=sp_data columns="area_sqm INT"
v.to.db sp_data option=area column=area_sqm unit=meters
Now convert to raster, using the area_sqm column as the raster value.
But first set the GRASS computational region:
g.region -p vector=sp_data res=100
Now create a vector grid of polygons at the raster resolution:
v.mkgrid map=grid
Intersect with the forest polygons, and calculate areas of these intersect polygons:
Intersect with the forest polygons, dissolve the polygons inside the same grid cell, calculate areas of these intersect polygons, and move the area value to the grid cells:
v.overlay ain=grid bin=sp_data operator=and output=forest_grid_intersect
v.dissolve --overwrite input=forest_grid_intersect@jn column=a_cat output=forest_grid_intersect_dissolve
v.db.addtable map=forest_grid_intersect_dissolve columns="area_sqm INT"
v.to.db forest_grid_intersect_dissolve option=area column=area_sqm unit=meters
v.db.join map=grid@jn column=cat other_table=forest_grid_intersect_dissolve other_column=cat
v.db.addcolumn map=forest_grid_intersect columns="area_sqm INT"
v.to.db forest_grid_intersect option=area column=area_sqm unit=meters
and convert to raster, using the area attribute as the raster value:
v.to.rast input=grid type=area output=sp_data_rast use=attr attribute_column=area_sqm
And you're done. Hope I have it right this time :-) .
Best Answer
Two things you could check:
First make sure that your resolution is set correctly. This you do with the "Edit Current Grass Region" button. You mentioned that the original csv file had points every 40 km. If you want to stay with that, then assuming your points are in a projected coordinate system you should set the resolution to 40,000. If the points are Lat/Lon then the cell size should be about 0.4 (4/10ths of a degree).
Next, I doubt that you want to use v.to.rast. That's for making a raster with individual pixel values matching the individual point values. (BTW, when you set use=val and value=1, you're telling grass to set all pixels to a value of 1. Not what you want...)
A better option might be to import the points directly into a raster using r.in.xyz. Of course, you'll need to create separate rasters for each month. The "Column of data values" will be one of the columns in the ascii file with the data for one of the months. And for the "Statistic" you'll want "mean".