| 1 | #!/usr/bin/env python
|
|---|
| 2 |
|
|---|
| 3 | """
|
|---|
| 4 | MODULE: v.habitat.dem
|
|---|
| 5 |
|
|---|
| 6 | AUTHOR(S): Helmut Kudrnovsky <alectoria AT gmx at>
|
|---|
| 7 |
|
|---|
| 8 | PURPOSE: DEM and solar derived characteristics of habitats
|
|---|
| 9 |
|
|---|
| 10 | COPYRIGHT: (C) 2014 by the GRASS Development Team
|
|---|
| 11 |
|
|---|
| 12 | This program is free software under the GNU General Public
|
|---|
| 13 | License (>=v2). Read the file COPYING that comes with GRASS
|
|---|
| 14 | for details.
|
|---|
| 15 | """
|
|---|
| 16 |
|
|---|
| 17 | #%module
|
|---|
| 18 | #% description: Calculates DEM derived characteristics of habitats.
|
|---|
| 19 | #% keyword: vector
|
|---|
| 20 | #% keyword: raster
|
|---|
| 21 | #% keyword: terrain
|
|---|
| 22 | #% keyword: statistics
|
|---|
| 23 | #% keyword: sun
|
|---|
| 24 | #% keyword: zonal statistics
|
|---|
| 25 | #%end
|
|---|
| 26 |
|
|---|
| 27 | #%option G_OPT_R_ELEV
|
|---|
| 28 | #% key: elevation
|
|---|
| 29 | #% description: Name of elevation raster map
|
|---|
| 30 | #% required: yes
|
|---|
| 31 | #%end
|
|---|
| 32 |
|
|---|
| 33 | #%option G_OPT_V_INPUT
|
|---|
| 34 | #% key: vector
|
|---|
| 35 | #% description: Name of habitat vector map
|
|---|
| 36 | #% required: yes
|
|---|
| 37 | #%end
|
|---|
| 38 |
|
|---|
| 39 | #%option G_OPT_DB_COLUMN
|
|---|
| 40 | #% description: Name of attribute column with a unique habitat ID (must be numeric)
|
|---|
| 41 | #% required: yes
|
|---|
| 42 | #%end
|
|---|
| 43 |
|
|---|
| 44 | #%option
|
|---|
| 45 | #% key: prefix
|
|---|
| 46 | #% type: string
|
|---|
| 47 | #% key_desc: prefix
|
|---|
| 48 | #% description: output prefix (must start with a letter)
|
|---|
| 49 | #% required: yes
|
|---|
| 50 | #%end
|
|---|
| 51 |
|
|---|
| 52 | #%option G_OPT_M_DIR
|
|---|
| 53 | #% key: dir
|
|---|
| 54 | #% description: Directory where the output will be found
|
|---|
| 55 | #% required : yes
|
|---|
| 56 | #%end
|
|---|
| 57 |
|
|---|
| 58 | #%option
|
|---|
| 59 | #% key: region_extension
|
|---|
| 60 | #% type: double
|
|---|
| 61 | #% key_desc: float
|
|---|
| 62 | #% description: region extension
|
|---|
| 63 | #% required : no
|
|---|
| 64 | #% answer: 5000
|
|---|
| 65 | #%end
|
|---|
| 66 |
|
|---|
| 67 | #%option
|
|---|
| 68 | #% key: start_time
|
|---|
| 69 | #% type: double
|
|---|
| 70 | #% label: Start time of interval
|
|---|
| 71 | #% description: Use up to 2 decimal places
|
|---|
| 72 | #% options: 0-24
|
|---|
| 73 | #% answer: 8
|
|---|
| 74 | #%end
|
|---|
| 75 |
|
|---|
| 76 | #%option
|
|---|
| 77 | #% key: end_time
|
|---|
| 78 | #% type: double
|
|---|
| 79 | #% label: End time of interval
|
|---|
| 80 | #% description: Use up to 2 decimal places
|
|---|
| 81 | #% options: 0-24
|
|---|
| 82 | #% answer: 18
|
|---|
| 83 | #%end
|
|---|
| 84 |
|
|---|
| 85 | #%option
|
|---|
| 86 | #% key: time_step
|
|---|
| 87 | #% type: double
|
|---|
| 88 | #% label: Time step for running r.sun [decimal hours]
|
|---|
| 89 | #% description: Use up to 2 decimal places
|
|---|
| 90 | #% options: 0-24
|
|---|
| 91 | #% answer: 1
|
|---|
| 92 | #%end
|
|---|
| 93 |
|
|---|
| 94 | #%option
|
|---|
| 95 | #% key: day
|
|---|
| 96 | #% type: integer
|
|---|
| 97 | #% description: No. of day of the year
|
|---|
| 98 | #% options: 1-365
|
|---|
| 99 | #% answer: 172
|
|---|
| 100 | #%end
|
|---|
| 101 |
|
|---|
| 102 | #%option
|
|---|
| 103 | #% key: year
|
|---|
| 104 | #% type: integer
|
|---|
| 105 | #% label: Year used for map registration into temporal dataset or r.timestamp
|
|---|
| 106 | #% description: This value is not used in r.sun calculations
|
|---|
| 107 | #% options: 1900-9999
|
|---|
| 108 | #% answer: 2014
|
|---|
| 109 | #%end
|
|---|
| 110 |
|
|---|
| 111 | import sys
|
|---|
| 112 | import os
|
|---|
| 113 | import grass.script as grass
|
|---|
| 114 | import math
|
|---|
| 115 | from numpy import array
|
|---|
| 116 | from numpy import zeros
|
|---|
| 117 | import csv
|
|---|
| 118 |
|
|---|
| 119 | def main():
|
|---|
| 120 | r_elevation = options['elevation'].split('@')[0]
|
|---|
| 121 | v_habitat = options['vector'].split('@')[0]
|
|---|
| 122 | v_column = options['column']
|
|---|
| 123 | regext = options['region_extension']
|
|---|
| 124 | directory = options['dir']
|
|---|
| 125 | prefix = options['prefix']
|
|---|
| 126 | r_accumulation = prefix+'_accumulation'
|
|---|
| 127 | r_drainage = prefix+'_drainage'
|
|---|
| 128 | r_tci = prefix+'_topographicindex'
|
|---|
| 129 | r_slope = prefix+'_slope'
|
|---|
| 130 | r_aspect = prefix+'_aspect'
|
|---|
| 131 | r_flow_accum = prefix+'_flow_accumulation'
|
|---|
| 132 | r_LS = prefix+'_LS_factor'
|
|---|
| 133 | r_habitat = prefix+'_'+v_habitat+'_rast'
|
|---|
| 134 | r_geomorphon = prefix+'_'+r_elevation+'_geomorphon'
|
|---|
| 135 | r_beam_rad = prefix+'_beam_rad'
|
|---|
| 136 | f_habitat_small_areas = prefix+'_small_areas.csv'
|
|---|
| 137 | f_habitat_geomorphons = prefix+'_habitat_geomorphons.csv'
|
|---|
| 138 | f_habitat_geomorphons_pivot = prefix+'_habitat_geomorphons_pivot.csv'
|
|---|
| 139 | f_LS_colorrules = prefix+'_LS_color.txt'
|
|---|
| 140 | t_habitat_geomorphons = prefix+'_habgeo_tab'
|
|---|
| 141 | t_habitat_geomorphons_pivot = prefix+'_habgeo_tab_pivot'
|
|---|
| 142 | d_start_time = options['start_time']
|
|---|
| 143 | d_end_time = options['end_time']
|
|---|
| 144 | d_time_step = options['time_step']
|
|---|
| 145 | d_day = options['day']
|
|---|
| 146 | d_year = options['year']
|
|---|
| 147 | saved_region = prefix+'_region_ori'
|
|---|
| 148 | current_region = grass.region()
|
|---|
| 149 | X = float(regext)
|
|---|
| 150 | N = current_region["n"]
|
|---|
| 151 | S = current_region["s"]
|
|---|
| 152 | E = current_region["e"]
|
|---|
| 153 | W = current_region["w"]
|
|---|
| 154 | Yres = current_region['nsres']
|
|---|
| 155 | Xres = current_region['ewres']
|
|---|
| 156 | PixelArea = Xres * Yres
|
|---|
| 157 | global tmp
|
|---|
| 158 |
|
|---|
| 159 | ## check if r.geomorphon addon is installed
|
|---|
| 160 | if not grass.find_program('r.geomorphon', '--help'):
|
|---|
| 161 | grass.fatal(_("The 'r.geomorphon' addon was not found, install it first:") +
|
|---|
| 162 | "\n" +
|
|---|
| 163 | "g.extension r.geomorphon")
|
|---|
| 164 |
|
|---|
| 165 | ## check if r.geomorphon addon is installed
|
|---|
| 166 | if not grass.find_program('r.sun.hourly', '--help'):
|
|---|
| 167 | grass.fatal(_("The 'r.sun.hourly' addon was not found, install it first:") +
|
|---|
| 168 | "\n" +
|
|---|
| 169 | "g.extension r.sun.hourly")
|
|---|
| 170 |
|
|---|
| 171 | # Region settings
|
|---|
| 172 | grass.message( "Current region will be saved and extended." )
|
|---|
| 173 | grass.message( "----" )
|
|---|
| 174 |
|
|---|
| 175 | # Print and save current region
|
|---|
| 176 | grass.message( "Current region:" )
|
|---|
| 177 | grass.message( "n, s, e, w" )
|
|---|
| 178 | grass.message( [current_region[key] for key in "nsew"] )
|
|---|
| 179 | grass.run_command('g.region', save = saved_region, overwrite = True)
|
|---|
| 180 | grass.message( "Current region saved." )
|
|---|
| 181 | grass.message( "----" )
|
|---|
| 182 |
|
|---|
| 183 | # does vector map exist in CURRENT mapset?
|
|---|
| 184 | mapset = grass.gisenv()['MAPSET']
|
|---|
| 185 | exists = bool(grass.find_file(v_habitat, element='vector', mapset=mapset)['file'])
|
|---|
| 186 | if not exists:
|
|---|
| 187 | grass.fatal(_("Vector map <%s> not found in current mapset") % v_habitat)
|
|---|
| 188 |
|
|---|
| 189 | # Align region to elevation raster and habitat vector
|
|---|
| 190 | grass.message( "Align region to elevation raster and habitat vector ..." )
|
|---|
| 191 | grass.run_command('g.region', flags = 'a',
|
|---|
| 192 | rast = r_elevation,
|
|---|
| 193 | vect = v_habitat,
|
|---|
| 194 | align = r_elevation)
|
|---|
| 195 | grass.message( "Alignment done." )
|
|---|
| 196 |
|
|---|
| 197 | aligned_region = grass.region()
|
|---|
| 198 | Naligned = aligned_region["n"]
|
|---|
| 199 | Saligned = aligned_region["s"]
|
|---|
| 200 | Ealigned = aligned_region["e"]
|
|---|
| 201 | Waligned = aligned_region["w"]
|
|---|
| 202 |
|
|---|
| 203 | grass.message( "Aligned region:" )
|
|---|
| 204 | grass.message( "n, s, e, w" )
|
|---|
| 205 | grass.message( [aligned_region[key] for key in "nsew"] )
|
|---|
| 206 | grass.message( "----" )
|
|---|
| 207 |
|
|---|
| 208 | # Extend region
|
|---|
| 209 | grass.message( "Extend region by" )
|
|---|
| 210 | grass.message( regext )
|
|---|
| 211 | grass.message( "in all directions" )
|
|---|
| 212 |
|
|---|
| 213 | grass.run_command('g.region', n = Naligned+X,
|
|---|
| 214 | s = Saligned-X,
|
|---|
| 215 | e = Ealigned+X,
|
|---|
| 216 | w = Waligned-X)
|
|---|
| 217 | grass.message( "Region extension done." )
|
|---|
| 218 |
|
|---|
| 219 | extended_region = grass.region()
|
|---|
| 220 |
|
|---|
| 221 | grass.message( "Extended region:" )
|
|---|
| 222 | grass.message( "n, s, e, w" )
|
|---|
| 223 | grass.message( [extended_region[key] for key in "nsew"] )
|
|---|
| 224 | grass.message( "----" )
|
|---|
| 225 |
|
|---|
| 226 | # Watershed calculation: accumulation, drainage direction, topographic index
|
|---|
| 227 | grass.message( "Calculation of accumulation, drainage direction, topographic index by r.watershed ..." )
|
|---|
| 228 | grass.run_command('r.watershed', elevation = r_elevation,
|
|---|
| 229 | accumulation = r_accumulation,
|
|---|
| 230 | drainage = r_drainage,
|
|---|
| 231 | tci = r_tci,
|
|---|
| 232 | convergence = 5,
|
|---|
| 233 | flags = 'am',
|
|---|
| 234 | overwrite = True)
|
|---|
| 235 | grass.message( "Calculation of accumulation, drainage direction, topographic done." )
|
|---|
| 236 | grass.message( "----" )
|
|---|
| 237 |
|
|---|
| 238 | # Calculation of slope and aspect maps
|
|---|
| 239 | grass.message( "Calculation of slope and aspect by r.slope.aspect ..." )
|
|---|
| 240 | grass.run_command('r.slope.aspect', elevation = r_elevation,
|
|---|
| 241 | slope = r_slope,
|
|---|
| 242 | aspect = r_aspect,
|
|---|
| 243 | overwrite = True)
|
|---|
| 244 | grass.message( "Calculation of slope and aspect done." )
|
|---|
| 245 | grass.message( "----" )
|
|---|
| 246 |
|
|---|
| 247 | # Calculate pixel area by nsres x ewres
|
|---|
| 248 | grass.message( "Pixel area:" )
|
|---|
| 249 | grass.message( PixelArea )
|
|---|
| 250 | grass.message( "----" )
|
|---|
| 251 |
|
|---|
| 252 | # Calculate habitat area and populate it to the attribute table
|
|---|
| 253 | grass.message( "Calculate habitat's areas and populate it to the attribute table ..." )
|
|---|
| 254 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 255 | layer = 1,
|
|---|
| 256 | columns = "habarea double")
|
|---|
| 257 |
|
|---|
| 258 | grass.run_command("v.to.db", map = v_habitat,
|
|---|
| 259 | option = 'area',
|
|---|
| 260 | layer = 1,
|
|---|
| 261 | columns = 'habarea',
|
|---|
| 262 | overwrite = True)
|
|---|
| 263 |
|
|---|
| 264 | grass.message( "Calculate habitat's areas done." )
|
|---|
| 265 | grass.message( "----" )
|
|---|
| 266 |
|
|---|
| 267 | # Show habitat areas smaller than Pixel Area
|
|---|
| 268 | grass.message( "Habitat areas smaller than pixel area." )
|
|---|
| 269 | grass.run_command("v.db.select", map = v_habitat,
|
|---|
| 270 | flags = 'v',
|
|---|
| 271 | layer = 1,
|
|---|
| 272 | columns = v_column,
|
|---|
| 273 | where = "habarea < %s" % (PixelArea))
|
|---|
| 274 |
|
|---|
| 275 | smallareacsv = os.path.join( directory, f_habitat_small_areas )
|
|---|
| 276 |
|
|---|
| 277 | grass.run_command("v.db.select", map = v_habitat,
|
|---|
| 278 | flags = 'v',
|
|---|
| 279 | layer = 1,
|
|---|
| 280 | columns = v_column,
|
|---|
| 281 | file = smallareacsv,
|
|---|
| 282 | where = "habarea < %s" % (PixelArea))
|
|---|
| 283 |
|
|---|
| 284 | grass.message( "A list of habitat areas smaller than pixel area can be found in: " )
|
|---|
| 285 | grass.message( smallareacsv )
|
|---|
| 286 | grass.message( "----" )
|
|---|
| 287 |
|
|---|
| 288 | # Mark habitats smaller than pixel area in attribute table
|
|---|
| 289 | grass.message( "Mark habitats smaller than pixel area in attribute table ..." )
|
|---|
| 290 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 291 | layer = 1,
|
|---|
| 292 | columns = "%s_smallarea varchar(1)" % prefix )
|
|---|
| 293 | grass.run_command("v.db.update", map = v_habitat, layer = 1, column = '%s_smallarea' % (prefix), value = '*', where = 'habarea < %s' % (PixelArea))
|
|---|
| 294 | grass.message( "See column" )
|
|---|
| 295 | grass.message( '%s_smallarea' % prefix )
|
|---|
| 296 | grass.message( "marked by *." )
|
|---|
| 297 | grass.message( "----" )
|
|---|
| 298 |
|
|---|
| 299 | # Upload DEM zonal statistics to the attribute table
|
|---|
| 300 | grass.message( "Upload DEM zonal statistics to the attribute table ..." )
|
|---|
| 301 | grass.run_command("v.rast.stats", map = v_habitat,
|
|---|
| 302 | flags = 'c',
|
|---|
| 303 | layer = 1,
|
|---|
| 304 | raster = r_elevation,
|
|---|
| 305 | column_prefix = prefix+'_dem',
|
|---|
| 306 | method = 'minimum,maximum,range,average,median')
|
|---|
| 307 | grass.message( "Upload DEM zonal statistics done." )
|
|---|
| 308 | grass.message( "----" )
|
|---|
| 309 |
|
|---|
| 310 | # Upload slope zonal statistics to the attribute table
|
|---|
| 311 | grass.message( "Upload slope zonal statistics to the attribute table ..." )
|
|---|
| 312 | grass.run_command("v.rast.stats", map = v_habitat,
|
|---|
| 313 | flags = 'c',
|
|---|
| 314 | layer = 1,
|
|---|
| 315 | raster = r_slope,
|
|---|
| 316 | column_prefix = prefix+'_slope',
|
|---|
| 317 | method = 'minimum,maximum,range,average,median')
|
|---|
| 318 | grass.message( "Upload slope zonal statistics done." )
|
|---|
| 319 | grass.message( "----" )
|
|---|
| 320 |
|
|---|
| 321 | # Upload slope zonal statistics to the attribute table
|
|---|
| 322 | grass.message( "Upload aspect zonal statistics to the attribute table ..." )
|
|---|
| 323 | grass.run_command("v.rast.stats", map = v_habitat,
|
|---|
| 324 | flags = 'c',
|
|---|
| 325 | layer = 1,
|
|---|
| 326 | raster = r_aspect,
|
|---|
| 327 | column_prefix = prefix+'_aspect',
|
|---|
| 328 | method = 'minimum,maximum,range,average,median')
|
|---|
| 329 | grass.message( "Upload aspect zonal statistics done." )
|
|---|
| 330 | grass.message( "----" )
|
|---|
| 331 |
|
|---|
| 332 | # Do some simple checks regarding aspect range
|
|---|
| 333 | grass.message( "Do some simple checks regarding aspect range and populate it to the attribute table..." )
|
|---|
| 334 | grass.message( "aspect range 100-200 *" )
|
|---|
| 335 | grass.message( "aspect range 201-300 **" )
|
|---|
| 336 | grass.message( "aspect range >= 300 ***" )
|
|---|
| 337 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 338 | layer = 1,
|
|---|
| 339 | columns = "%s varchar(3)" % (prefix+'_check_aspect_range'))
|
|---|
| 340 |
|
|---|
| 341 | grass.run_command("db.execute", sql = "UPDATE %s SET %s ='*' WHERE %s < 200 AND %s >= 100" % (v_habitat, prefix+'_check_aspect_range', prefix+'_aspect_range', prefix+'_aspect_range'))
|
|---|
| 342 | grass.run_command("db.execute", sql = "UPDATE %s SET %s ='**' WHERE %s < 300 AND %s >= 200" % (v_habitat, prefix+'_check_aspect_range', prefix+'_aspect_range', prefix+'_aspect_range'))
|
|---|
| 343 | grass.run_command("db.execute", sql = "UPDATE %s SET %s ='***' WHERE %s >= 300" % (v_habitat, prefix+'_check_aspect_range', prefix+'_aspect_range'))
|
|---|
| 344 |
|
|---|
| 345 |
|
|---|
| 346 | grass.message( "Simple checks regarding aspect range done." )
|
|---|
| 347 | grass.message( "----" )
|
|---|
| 348 |
|
|---|
| 349 | # Do some simple checks regarding aspect and and slope
|
|---|
| 350 | grass.message( "Do some simple checks regarding aspect range and slope median and populate it to the attribute table..." )
|
|---|
| 351 | grass.message( "aspect range 100-200 and median slope < 5 *" )
|
|---|
| 352 | grass.message( "aspect range 201-300 and median slope < 5 **" )
|
|---|
| 353 | grass.message( "aspect range >= 300 and median slope < 5 ***" )
|
|---|
| 354 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 355 | layer = 1,
|
|---|
| 356 | columns = "%s varchar(3)" % (prefix+'_check_aspect_slope'))
|
|---|
| 357 |
|
|---|
| 358 | grass.run_command("db.execute", sql = "UPDATE %s SET %s ='*' WHERE (%s < 200 AND %s >= 100) AND %s < 5" % (v_habitat, prefix+'_check_aspect_slope', prefix+'_aspect_range', prefix+'_aspect_range', prefix+'_slope_median'))
|
|---|
| 359 | grass.run_command("db.execute", sql = "UPDATE %s SET %s ='**' WHERE (%s < 300 AND %s >= 200) AND %s < 5" % (v_habitat, prefix+'_check_aspect_slope', prefix+'_aspect_range', prefix+'_aspect_range', prefix+'_slope_median'))
|
|---|
| 360 | grass.run_command("db.execute", sql = "UPDATE %s SET %s ='***' WHERE %s >= 300 AND %s < 5" % (v_habitat, prefix+'_check_aspect_slope', prefix+'_aspect_range', prefix+'_slope_median'))
|
|---|
| 361 |
|
|---|
| 362 | grass.message( "Simple checks regarding aspect range and median slope done." )
|
|---|
| 363 | grass.message( "----" )
|
|---|
| 364 |
|
|---|
| 365 | # Do some simple checks regarding aspect and and slope
|
|---|
| 366 | grass.message( "Convert habitat vector to raster ..." )
|
|---|
| 367 | grass.run_command("v.to.rast", input = v_habitat,
|
|---|
| 368 | layer = 1,
|
|---|
| 369 | type = 'area',
|
|---|
| 370 | use = 'attr',
|
|---|
| 371 | attrcolumn = v_column,
|
|---|
| 372 | output = r_habitat)
|
|---|
| 373 | grass.message( "Conversion done." )
|
|---|
| 374 | grass.message( "----" )
|
|---|
| 375 |
|
|---|
| 376 | # Calculate the most common geomorphons
|
|---|
| 377 | grass.message( "Calculate the most common geomorphons ..." )
|
|---|
| 378 | grass.run_command("r.geomorphon", elevation = r_elevation,
|
|---|
| 379 | skip = 0,
|
|---|
| 380 | search = 3,
|
|---|
| 381 | flat = 1,
|
|---|
| 382 | dist = 0,
|
|---|
| 383 | forms = r_geomorphon)
|
|---|
| 384 | grass.message( "Geomorphon calculations done." )
|
|---|
| 385 | grass.message( "----" )
|
|---|
| 386 |
|
|---|
| 387 | # Mutual occurrence (coincidence) of categories of habitats and geomorphons
|
|---|
| 388 | grass.message( "Calculate mutual occurrences of habitats and geomorphons ..." )
|
|---|
| 389 | grass.message( "1 - flat" )
|
|---|
| 390 | grass.message( "2 - summit" )
|
|---|
| 391 | grass.message( "3 - ridge" )
|
|---|
| 392 | grass.message( "4 - shoulder" )
|
|---|
| 393 | grass.message( "5 - spur" )
|
|---|
| 394 | grass.message( "6 - slope" )
|
|---|
| 395 | grass.message( "7 - hollow" )
|
|---|
| 396 | grass.message( "8 - footslope" )
|
|---|
| 397 | grass.message( "9 - valley" )
|
|---|
| 398 | grass.message( "10 - depression" )
|
|---|
| 399 | grass.message( " " )
|
|---|
| 400 | grass.message( "Mutual occurrence in percent of the row" )
|
|---|
| 401 | grass.run_command("r.coin", first = r_habitat,
|
|---|
| 402 | second = r_geomorphon,
|
|---|
| 403 | flags = 'w',
|
|---|
| 404 | units = 'y')
|
|---|
| 405 | grass.message( "Calculations of mutual occurrences done." )
|
|---|
| 406 | grass.message( "----" )
|
|---|
| 407 |
|
|---|
| 408 | # Join geomorphons to habitat attribute table
|
|---|
| 409 | grass.message( "Join geomorphon information to habitat attribute table ...." )
|
|---|
| 410 |
|
|---|
| 411 | habgeocsv = os.path.join( directory, f_habitat_geomorphons)
|
|---|
| 412 |
|
|---|
| 413 | grass.run_command("r.stats", input = [r_habitat, r_geomorphon],
|
|---|
| 414 | flags = 'aln',
|
|---|
| 415 | separator = ';',
|
|---|
| 416 | output = habgeocsv)
|
|---|
| 417 |
|
|---|
| 418 | grass.run_command("db.in.ogr", input = habgeocsv,
|
|---|
| 419 | output = t_habitat_geomorphons)
|
|---|
| 420 |
|
|---|
| 421 | grass.run_command("db.dropcolumn", table = t_habitat_geomorphons,
|
|---|
| 422 | column = 'field_2',
|
|---|
| 423 | flags = 'f')
|
|---|
| 424 |
|
|---|
| 425 | grass.run_command("db.dropcolumn", table = t_habitat_geomorphons,
|
|---|
| 426 | column = 'field_3',
|
|---|
| 427 | flags = 'f')
|
|---|
| 428 |
|
|---|
| 429 | habgeocsv_pivot = os.path.join( directory, f_habitat_geomorphons_pivot)
|
|---|
| 430 |
|
|---|
| 431 | grass.run_command("db.select", separator = ';',
|
|---|
| 432 | output = habgeocsv_pivot,
|
|---|
| 433 | sql = "SELECT field_1, sum(case when field_4 = 'flat' then field_5 end) as flat, sum(case when field_4 = 'summit' then field_5 end) as summit, sum(case when field_4 = 'ridge' then field_5 end) as ridge, sum(case when field_4 = 'shoulder' then field_5 end) as shoulder, sum(case when field_4 = 'spur' then field_5 end) as spur, sum(case when field_4 = 'slope' then field_5 end) as slope, sum(case when field_4 = 'hollow' then field_5 end) as hollow, sum(case when field_4 = 'footslope' then field_5 end) as footslope, sum(case when field_4 = 'valley' then field_5 end) as valley, sum(case when field_4 = 'depression' then field_5 end) as depression , sum(field_5) as SubTotal FROM %s GROUP BY field_1" % t_habitat_geomorphons)
|
|---|
| 434 |
|
|---|
| 435 | grass.run_command("db.in.ogr", input = habgeocsv_pivot,
|
|---|
| 436 | output = t_habitat_geomorphons_pivot)
|
|---|
| 437 |
|
|---|
| 438 | grass.run_command("v.db.join", map = v_habitat,
|
|---|
| 439 | column = v_column,
|
|---|
| 440 | other_table = t_habitat_geomorphons_pivot,
|
|---|
| 441 | other_column = 'field_1')
|
|---|
| 442 | grass.message( "Geomorphon information joint to habitat attribute table." )
|
|---|
| 443 | grass.message( "..." )
|
|---|
| 444 | grass.message( "Calculating percent geomorphon of habitat area ..." )
|
|---|
| 445 | # add column for percent geomorphon
|
|---|
| 446 |
|
|---|
| 447 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 448 | layer = 1,
|
|---|
| 449 | columns = "%s_perc_flat double precision" % prefix )
|
|---|
| 450 |
|
|---|
| 451 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 452 | layer = 1,
|
|---|
| 453 | columns = "%s_perc_summit double precision" % prefix )
|
|---|
| 454 |
|
|---|
| 455 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 456 | layer = 1,
|
|---|
| 457 | columns = "%s_perc_ridge double precision" % prefix )
|
|---|
| 458 |
|
|---|
| 459 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 460 | layer = 1,
|
|---|
| 461 | columns = "%s_perc_shoulder double precision" % prefix )
|
|---|
| 462 |
|
|---|
| 463 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 464 | layer = 1,
|
|---|
| 465 | columns = "%s_perc_spur double precision" % prefix )
|
|---|
| 466 |
|
|---|
| 467 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 468 | layer = 1,
|
|---|
| 469 | columns = "%s_perc_slope double precision" % prefix )
|
|---|
| 470 |
|
|---|
| 471 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 472 | layer = 1,
|
|---|
| 473 | columns = "%s_perc_hollow double precision" % prefix )
|
|---|
| 474 |
|
|---|
| 475 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 476 | layer = 1,
|
|---|
| 477 | columns = "%s_perc_footslope double precision" % prefix )
|
|---|
| 478 |
|
|---|
| 479 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 480 | layer = 1,
|
|---|
| 481 | columns = "%s_perc_valley double precision" % prefix )
|
|---|
| 482 |
|
|---|
| 483 | grass.run_command("v.db.addcolumn", map = v_habitat,
|
|---|
| 484 | layer = 1,
|
|---|
| 485 | columns = "%s_perc_depression double precision" % prefix )
|
|---|
| 486 |
|
|---|
| 487 | # calculate percent geomorphon
|
|---|
| 488 |
|
|---|
| 489 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 490 | layer = 1,
|
|---|
| 491 | column = "%s_perc_flat" % prefix,
|
|---|
| 492 | query_column="cast(flat AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 493 |
|
|---|
| 494 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 495 | layer = 1,
|
|---|
| 496 | column = "%s_perc_summit" % prefix,
|
|---|
| 497 | query_column="cast(summit AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 498 |
|
|---|
| 499 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 500 | layer = 1,
|
|---|
| 501 | column = "%s_perc_ridge" % prefix,
|
|---|
| 502 | query_column="cast(ridge AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 503 |
|
|---|
| 504 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 505 | layer = 1,
|
|---|
| 506 | column = "%s_perc_shoulder" % prefix,
|
|---|
| 507 | query_column="cast(shoulder AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 508 |
|
|---|
| 509 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 510 | layer = 1,
|
|---|
| 511 | column = "%s_perc_spur" % prefix,
|
|---|
| 512 | query_column="cast(spur AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 513 |
|
|---|
| 514 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 515 | layer = 1,
|
|---|
| 516 | column = "%s_perc_slope" % prefix,
|
|---|
| 517 | query_column="cast(slope AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 518 |
|
|---|
| 519 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 520 | layer = 1,
|
|---|
| 521 | column = "%s_perc_hollow" % prefix,
|
|---|
| 522 | query_column="cast(hollow AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 523 |
|
|---|
| 524 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 525 | layer = 1,
|
|---|
| 526 | column = "%s_perc_footslope" % prefix,
|
|---|
| 527 | query_column="cast(footslope AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 528 |
|
|---|
| 529 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 530 | layer = 1,
|
|---|
| 531 | column = "%s_perc_valley" % prefix,
|
|---|
| 532 | query_column="cast(valley AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 533 |
|
|---|
| 534 | grass.run_command("v.db.update", map = v_habitat,
|
|---|
| 535 | layer = 1,
|
|---|
| 536 | column = "%s_perc_depression" % prefix,
|
|---|
| 537 | query_column="cast(depression AS real) / cast( SubTotal AS real) * 100.0")
|
|---|
| 538 |
|
|---|
| 539 | grass.message( " " )
|
|---|
| 540 | grass.message( "Calculating of percent geomorphon of habitat area done." )
|
|---|
| 541 | grass.message( "----" )
|
|---|
| 542 |
|
|---|
| 543 |
|
|---|
| 544 | # Give information where output files are
|
|---|
| 545 | grass.message( "Geomorphon information:" )
|
|---|
| 546 | grass.message( habgeocsv )
|
|---|
| 547 | grass.message( "Geomorphon information in pivot format:" )
|
|---|
| 548 | grass.message( habgeocsv_pivot )
|
|---|
| 549 | grass.message( "----" )
|
|---|
| 550 |
|
|---|
| 551 | # Calculate LS factor see Neteler & Mitasova 2008. Open Source GIS - A GRASS GIS Approach
|
|---|
| 552 | grass.message( "Calculate LS factor ..." )
|
|---|
| 553 | grass.run_command("r.flow", elevation = r_elevation,
|
|---|
| 554 | aspect = r_aspect,
|
|---|
| 555 | flowaccumulation = r_flow_accum)
|
|---|
| 556 |
|
|---|
| 557 | grass.message( "..." )
|
|---|
| 558 | grass.mapcalc("$outmap = 1.4 * exp($flowacc * $resolution / 22.1, 0.4) * exp(sin($slope) / 0.09, 1.2)",
|
|---|
| 559 | outmap = r_LS,
|
|---|
| 560 | flowacc = r_flow_accum,
|
|---|
| 561 | resolution = Xres,
|
|---|
| 562 | slope = r_slope)
|
|---|
| 563 |
|
|---|
| 564 | # create and define color rules file for LS factor map
|
|---|
| 565 | ls_color_rules_out = os.path.join( directory, f_LS_colorrules )
|
|---|
| 566 | with open(ls_color_rules_out, 'w') as f:
|
|---|
| 567 | writer = csv.writer(f)
|
|---|
| 568 | writer.writerow(['0 white'])
|
|---|
| 569 | writer.writerow(['3 yellow'])
|
|---|
| 570 | writer.writerow(['6 orange'])
|
|---|
| 571 | writer.writerow(['10 red'])
|
|---|
| 572 | writer.writerow(['50 magenta'])
|
|---|
| 573 | writer.writerow(['100 violet'])
|
|---|
| 574 |
|
|---|
| 575 | grass.run_command("r.colors", map = r_LS,
|
|---|
| 576 | rules = ls_color_rules_out)
|
|---|
| 577 |
|
|---|
| 578 | grass.message( "Calculation LS factor done." )
|
|---|
| 579 | grass.message( "----" )
|
|---|
| 580 |
|
|---|
| 581 | # Run r.sun.hourly in binary mode for light/shadow
|
|---|
| 582 | grass.message( "Run r.sun.hourly in binary mode for light/shadow for a certain day in the year ..." )
|
|---|
| 583 | grass.run_command("r.sun.hourly", elevation = r_elevation,
|
|---|
| 584 | flags = 'tb',
|
|---|
| 585 | aspect = r_aspect,
|
|---|
| 586 | slope = r_slope,
|
|---|
| 587 | start_time = d_start_time,
|
|---|
| 588 | end_time = d_end_time,
|
|---|
| 589 | day = d_day,
|
|---|
| 590 | year = d_year,
|
|---|
| 591 | beam_rad_basename = r_beam_rad)
|
|---|
| 592 |
|
|---|
| 593 | grass.message( "----" )
|
|---|
| 594 | grass.message( "Light/shadow conditions calculated for year" )
|
|---|
| 595 | grass.message( d_year )
|
|---|
| 596 | grass.message( "and day" )
|
|---|
| 597 | grass.message( d_day )
|
|---|
| 598 | grass.message( 'from' )
|
|---|
| 599 | grass.message( d_start_time )
|
|---|
| 600 | grass.message( 'to' )
|
|---|
| 601 | grass.message( d_end_time )
|
|---|
| 602 | grass.message( 'done.' )
|
|---|
| 603 | grass.message( "----" )
|
|---|
| 604 | grass.run_command("t.info", flags = 'h',
|
|---|
| 605 | input = r_beam_rad)
|
|---|
| 606 | grass.message( "----" )
|
|---|
| 607 |
|
|---|
| 608 | # Set region to original
|
|---|
| 609 | grass.message( "Restore original region settings:" )
|
|---|
| 610 | grass.run_command("g.region", flags = 'p',
|
|---|
| 611 | region = saved_region)
|
|---|
| 612 | grass.message( "----" )
|
|---|
| 613 |
|
|---|
| 614 | # clean up some temporay files and maps
|
|---|
| 615 | grass.message( "Some clean up ..." )
|
|---|
| 616 | grass.run_command("g.remove", flags="f", type="region", name= saved_region)
|
|---|
| 617 | grass.run_command("g.remove", flags="f", type="raster", name= r_flow_accum)
|
|---|
| 618 | grass.run_command("g.remove", flags="f", type="raster", name= r_habitat)
|
|---|
| 619 | grass.run_command("db.droptable", flags = 'f',
|
|---|
| 620 | table = t_habitat_geomorphons)
|
|---|
| 621 | grass.run_command("db.droptable", flags = 'f',
|
|---|
| 622 | table = t_habitat_geomorphons_pivot)
|
|---|
| 623 | grass.run_command("db.dropcolumn", flags = 'f',
|
|---|
| 624 | table = v_habitat,
|
|---|
| 625 | column = 'field_1')
|
|---|
| 626 | grass.message( "Clean up done." )
|
|---|
| 627 | grass.message( "----" )
|
|---|
| 628 |
|
|---|
| 629 | # v.habitat.dem done!
|
|---|
| 630 | grass.message( "v.habitat.dem done!" )
|
|---|
| 631 |
|
|---|
| 632 | if __name__ == "__main__":
|
|---|
| 633 | options, flags = grass.parser()
|
|---|
| 634 | sys.exit(main())
|
|---|
