In GRASS GIS:
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.
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.Sign up to join this community
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:
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