I have a shapefile with 300 hundred of sites (with height values for each point), and want to calculate the percentage or number of cells (of a given DEM) which are lower than the site's value.

  • Can the answer use bash or python scripting? Also, do the data need to be output in a particular way, e.g., to text file or updated as an attribute to the original shapefile? – user55937 Feb 2 '16 at 19:24
  • Actually I don't know python scripting, but someone can help in this. It would be fine if the output is an attribute to the original file. I hoped the answer could be doing something with the raster calculator... – user66349 Feb 2 '16 at 19:37
  • I posted some code. Let me know if this doesn't satisfy your requirements. Also, you said you had some to help you get up and running with python, but let me know if I can help with that too. – user55937 Feb 2 '16 at 23:49
  • @user55937 , the person who I supposed to help is unable to do that. Can you help me to run this code on GRASS? – user66349 Feb 4 '16 at 9:06
  • 1
    The easiest way to use the script is probably to call python from inside the grass shell, and then copy and paste the code into the interpreter. That is, when you start grass, you usually have some type of shell also. Type 'python' into the shell and hit enter. You should see python start up. If this doesn't work, let me know. And if not, let me know what os you're using. – user55937 Feb 4 '16 at 11:11

Here's a small python script. The DEM is from the GRASS lidar time series dataset (https://grass.osgeo.org/download/sample-data/). The contents of pts.txt referenced in the script (I assume this is like your sites point data) are:


Here's the points on the DEM: enter image description here Here's the script that uses mapcalc to calculate the percentage of DEM under the elevation of the site locations:

import grass.script as grass
from string import digits

points_file = 'pts.txt'
dem_name = 'JR_NH_20080327_1m'

grass.run_command('v.in.ascii', input=points_file, output='sites', format='point', fs=',', x=1, y=2)

site_locations = grass.read_command('v.out.ascii', input='sites', output='-', format='point', fs=',')
site_locations = site_locations.strip('\n').split('\n')

for site_location in site_locations:
    # remove category from site_location: 913606.791589,250337.104673,1 to 913606.791589,250337.104673
    site_location = site_location.rstrip(digits)
    site_location = site_location.rstrip(',')
    # query raster with r.what
    site_query = grass.read_command('r.what', input=dem_name, east_north=site_location) 
    # get z value from r.what result: '913606.791589|250337.104673||5.937129\n'
    site_elevation = float( site_query.split('|')[-1] )
    # use mapcalc to get a map where 1 is under site and zero is above site
    exp = 'temp=' + dem_name + '<' + str(site_elevation)
    grass.mapcalc(exp, overwrite=True)
    # use r.univar to get mean, which the is percent lower than site
    rast_stats = grass.parse_command('r.univar',flags='g',map='temp')
    m = float(rast_stats['mean'])
    print str(m*100) + '% of dem is under ' + str(site_elevation)

which produces the output:

66.4153955206% of dem is under 5.937129
6.616916404% of dem is under 1.394099

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.