Commit 29cf249d authored by Roelofsen, Hans's avatar Roelofsen, Hans
Browse files

MNP experimenting scripts

parent 3e3836cd
###############################################################################
# Generate Habitat Suitability raster for pH for a single species
import rasterio as rio
import numpy as np
# Read ph raster from file
ph = rio.open(r'c:\apps\temp_geodata\zuurtegraad\ph_huidig_1990-2010_2016.flt')
# Get species response to pH.
# W:\PROJECTS\QMAR\MNP-SNL-ParameterSet\Parameters_v05_2019_06_12\06_MNP_versie4_par_response_PH.csv
# Species_code,Response_L20,Response_L80,Response_H80,Response_H20
# S02004210,4,5,7.17,9999
response_breakpoints = np.array([4, 5, 7.17, 9999])
response_capacity = [0, 0.5, 1, 0.5, 0]
mydict = dict(enumerate(response_capacity))
# Use numpy digitize: https://numpy.org/doc/stable/reference/generated/numpy.digitize.html
ph_response_array0 = np.vectorize(mydict.get)(np.digitize(ph.read(1), response_breakpoints))
###############################################################################
# Identify connected habitat areas and key areas
import numpy as np
import rasterio as rio
from skimage import measure
from skimage import morphology
# read dummy habitat raster. Must be 0/1 raster
habitat = rio.open(r'c:\apps\temp_geodata\sampleBT\N1001.tif')
habitat_arr = habitat.read(1)
assert habitat_arr.max(), habitat_arr.min() == (1, 0)
# Species local distance and key areas
species_key_area = 10 # hectare
species_local_distance_m = 100 # meter
# Get pixel size and translate local distance to pixels
pxl_size = int(habitat.res[0])
assert species_local_distance_m % pxl_size == 0
species_local_distance_pxls = int(species_local_distance_m / pxl_size)
# Generate search window to buffer habitats, with radius species_local_distance / 2
footprint = morphology.disk(int(species_local_distance_pxls / 2))
# Dilate habitat raster with species local distance: https://homepages.inf.ed.ac.uk/rbf/HIPR2/dilate.htm
dilated = morphology.dilation(habitat_arr, footprint)
# Calculated connected components based on dilated image. Use horizontal/vertical connectivity only
# https://scikit-image.org/docs/stable/api/skimage.measure.html#label
regions, nregions = measure.label(dilated, return_num=True, connectivity=1)
# Restrict labeled components to original habitat locations
habitat_regions = np.where(habitat_arr == 1, regions, 0)
# Get area of habitat regions and identify key areas.
vals, counts = np.unique(habitat_regions, return_counts=True)
habitat_areas_ha = zip(vals, np.multiply(counts, np.divide(np.square(pxl_size), 10000)))
key_areas = [k for k, v in habitat_areas_ha if v >= species_key_area and k > 0]
# Output profile. Find datatype to accomodate nregions
prof = habitat.profile
prof['dtype'] = rio.dtypes.get_minimum_dtype([0, nregions])
prof['nodata'] = 0
# Write connected areas to file.
with rio.open(r'c:\apps\temp_geodata\sampleBT\N1001_{}m_connected.tif'.format(species_local_distance_m), 'w', **prof) as dest:
dest.write(habitat_regions, 1)
# Write Key Areas to file.
with rio.open(r'c:\apps\temp_geodata\sampleBT\N1001_{}ha_KeyAreas.tif'.format(species_key_area), 'w', **prof) as dest:
dest.write(np.where(np.isin(habitat_regions, key_areas), 1, 0), 1)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment