Talk:FUTURES land-change modeling for evaluating innovative conservation scenarios: Difference between revisions
Jump to navigation
Jump to search
(Created page with "= Python: Loss of farmland and forest = <source lang="python"> import numpy as np import grass.script as gscript # list maps computed by r.futures.parallelpga for one scena...") |
⚠️Gmsanchez (talk | contribs) No edit summary |
||
(4 intermediate revisions by the same user not shown) | |||
Line 7: | Line 7: | ||
current_mapset = gscript.gisenv()['MAPSET'] | current_mapset = gscript.gisenv()['MAPSET'] | ||
# put the name you are using for the scenarios in the pattern option | # put the name you are using for the scenarios in the pattern option | ||
maps = gscript.list_grouped('raster', pattern=" | maps = gscript.list_grouped('raster', pattern="final*")[current_mapset] | ||
developed = [] | |||
forest = [] | forest = [] | ||
agriculture = [] | agriculture = [] | ||
Line 13: | Line 14: | ||
for each in maps: | for each in maps: | ||
# compute landcover in 2035 by updating landcover 2011 with the simulation results: | # compute landcover in 2035 by updating landcover 2011 with the simulation results: | ||
gscript.mapcalc("landuse_2035_{name} = if ({name} >=1, 21, landuse_2011)".format(name=each), quiet=True) | gscript.mapcalc("landuse_2035_{name} = if({name} >=1, 21, landuse_2011)".format(name=each), quiet=True, overwrite=True) | ||
# get information about how many cells there are of each land use category, use r.stats | # get information about how many cells there are of each land use category, use r.stats | ||
res = gscript.read_command('r.stats', flags='cn', input="landuse_2035_{name}".format(name=each), quiet=True).strip() | res = gscript.read_command('r.stats', flags='cn', input="landuse_2035_{name}".format(name=each), quiet=True).strip() | ||
Line 46: | Line 47: | ||
= Python: Forest fragmentation = | = Python: Forest fragmentation = | ||
We will need the addon r.forestfrag. You can install it in the command console with: | |||
<source lang="bash"> | |||
g.extension r.forestfrag | |||
</source> | |||
<source lang="python"> | <source lang="python"> | ||
import grass.script as gscript | import grass.script as gscript |
Latest revision as of 22:38, 3 April 2019
Python: Loss of farmland and forest
import numpy as np
import grass.script as gscript
# list maps computed by r.futures.parallelpga for one scenario
current_mapset = gscript.gisenv()['MAPSET']
# put the name you are using for the scenarios in the pattern option
maps = gscript.list_grouped('raster', pattern="final*")[current_mapset]
developed = []
forest = []
agriculture = []
other = []
for each in maps:
# compute landcover in 2035 by updating landcover 2011 with the simulation results:
gscript.mapcalc("landuse_2035_{name} = if({name} >=1, 21, landuse_2011)".format(name=each), quiet=True, overwrite=True)
# get information about how many cells there are of each land use category, use r.stats
res = gscript.read_command('r.stats', flags='cn', input="landuse_2035_{name}".format(name=each), quiet=True).strip()
for line in res.strip().splitlines():
cat, cells = line.split()
cat = int(cat)
cells = int(cells)
# developed
if 21 <= cat <= 24:
developed.append(cells)
# agriculture
elif cat == 81 or cat == 82:
agriculture.append(cells)
# forest
elif 41 <= cat <= 43:
forest.append(cells)
# other land use
else:
other.append(cells)
# convert to NumPy arrays
forest = np.array(forest)
agriculture = np.array(agriculture)
other = np.array(other)
developed = np.array(developed)
print "forest: {mean} +- {std}".format(mean=np.mean(forest), std=np.std(forest))
print "agriculture: {mean} +- {std}".format(mean=np.mean(agriculture), std=np.std(agriculture))
print "other: {mean} +- {std}".format(mean=np.mean(other), std=np.std(other))
print "developed: {mean} +- {std}".format(mean=np.mean(developed), std=np.std(developed))
Python: Forest fragmentation
We will need the addon r.forestfrag. You can install it in the command console with:
g.extension r.forestfrag
import grass.script as gscript
def forest_fragmentation(landuse):
gscript.mapcalc("forest_2035 = if({landuse} >= 41 && {landuse} <= 43, 1, 0)".format(landuse=landuse))
gscript.run_command('r.forestfrag', flags='a', input='forest_2035', output='fragment', window=25)
result = gscript.parse_command('r.stats', input='fragment', flags='cn', parse=(gscript.parse_key_val, {'sep': ' ', 'val_type': int}))
gscript.run_command('g.remove', type='raster', name=['forest_2035', 'fragment'], flags='f')
return result
current_mapset = gscript.gisenv()['MAPSET']
landuse_maps = gscript.list_grouped('raster', pattern="landuse_2035_*")[current_mapset]
landuse_maps = landuse_maps[:5]
for landuse in landuse_maps:
result = forest_fragmentation(landuse)
print result