This notebook creates a map of biophony and anthrophony, predicted from a multilevel regression model.
The model takes land cover areas (within a specified radius) as input, which this notebook computes for a grid of points over the study area. The output is then overlayed on web background tiles over the study area.
In [1]:
# datawaves database
from landscape.models import LandCoverMergedMapArea
from database.models import Site
from geo.models import Boundary
from django.contrib.gis.geos import Point, Polygon
from django.contrib.gis.db.models.functions import Intersection, Envelope
from django.contrib.gis.db.models import Collect
from os import path
import numpy
import pandas
import pyprind
import pyproj
import folium
from folium.plugins import ImageOverlay
from PIL import Image
from matplotlib import pyplot
from matplotlib import cm
%matplotlib inline
# hd5 to save results between sessions
import h5py
radius — the radius in meters for which to compute land cover areas around each grid location
spacing — the spacing in meters between each grid location
p1, p2 — projections defining the local coordinate system (p1) and the web mapping coordinate system (p2)
In [2]:
radius = 500
In [3]:
spacing = 200
In [4]:
p1 = pyproj.Proj(init='EPSG:31254') # MGI / Austria GK West
p2 = pyproj.Proj(init='EPSG:4326') # WGS 84
In [5]:
# naturalness = (%natural - %anthropogenic) / (% natural + % anthropogenic)
def get_naturalness(point, radius=radius):
buffer = point.buffer(radius)
result = LandCoverMergedMapArea.objects.filter(geometry__intersects=buffer, cover_type__in=[1, 2, 9])\
.annotate(intersection=Intersection('geometry', buffer))
forest = result.filter(cover_type__exact=9)
forest = forest.aggregate(total=Collect('intersection'))
result = result.aggregate(total=Collect('intersection'))
try:
forest_area = forest['total'].area
except AttributeError:
forest_area = 0
try:
result_area = result['total'].area
except AttributeError:
result_area = 0
urban_area = result_area - forest_area
try:
naturalness = (forest_area - urban_area) / (result_area)
except ZeroDivisionError:
naturalness = 0
return naturalness
In [6]:
# efficently queries the land cover data
# to determine the percentage of the specified land cover
# within a specified radius around a point location
def get_forest_net_area(point, radius=radius):
buffer = point.buffer(radius)
result = LandCoverMergedMapArea.objects.filter(shape__intersects=buffer, cover_type__in=[1, 2, 9])\
.annotate(intersection=Intersection('shape', buffer))
forest = result.filter(cover_type__exact=9)
forest = forest.aggregate(total=Collect('intersection'))
result = result.aggregate(total=Collect('intersection'))
try:
forest_area = forest['total'].area
except AttributeError:
forest_area = 0
try:
result_area = result['total'].area
except AttributeError:
result_area = 0
net_area = (2 * forest_area) - result_area
return (net_area / buffer.area) * 100
In [7]:
def get_pavement_area(point, radius=radius):
buffer = point.buffer(radius)
result = LandCoverMergedMapArea.objects.filter(geometry__intersects=buffer, cover_type__in=[1, 2, 9])\
.annotate(intersection=Intersection('geometry', buffer))
pavement = result.filter(cover_type__exact=2)
pavement = pavement.aggregate(total=Collect('intersection'))
try:
pavement_area = pavement['total'].area
except AttributeError:
pavement_area = 0
return (pavement_area / buffer.area) * 100
In [8]:
# utility function to save results between sessions
def export_to_hdf5(filepath, name, data):
h5f = h5py.File(filepath, 'w')
h5f.create_dataset(name, data=data)
h5f.close()
load study area boundary from the database
In [9]:
boundary = Boundary.objects.get(name = "study area").geometry.envelope
create a grid of point coordinates (in a numpy array) based on the boundary and desired node spacing
In [10]:
coords = numpy.array(boundary.coords[0])
bbox = dict()
bbox['left'] = coords[:, 0].min()
bbox['bottom'] = coords[:, 1].min()
bbox['right'] = coords[:, 0].max()
bbox['top'] = coords[:, 1].max()
rows = int(numpy.ceil((bbox['top'] - bbox['bottom']) / spacing))
columns = int(numpy.ceil((bbox['right'] - bbox['left']) / spacing))
grid = numpy.empty(shape=(rows, columns), dtype=numpy.ndarray)
x_start = bbox['left']
y_start = bbox['bottom']
for row in range(rows):
for column in range(columns):
grid[row, column] = [x_start + (spacing * column), y_start + (spacing * row)]
In [11]:
extended_bbox = dict()
extended_bbox['left'] = 71753.8726
extended_bbox['bottom'] = 229581.4586
extended_bbox['right'] = 86753.8726
extended_bbox['top'] = 241581.4586
extended_rows = int((extended_bbox['top'] - extended_bbox['bottom']) / spacing)
extended_columns = int((extended_bbox['right'] - extended_bbox['left']) / spacing)
extended_grid = numpy.empty(shape=(extended_rows, extended_columns), dtype=numpy.ndarray)
extended_x_start = extended_bbox['left']
extended_y_start = extended_bbox['bottom']
for row in range(extended_rows):
for column in range(extended_columns):
extended_grid[row, column] = [extended_x_start + (spacing * column), extended_y_start + (spacing * row)]
In [12]:
grid_extent = Polygon(((bbox['left'], bbox['bottom']),
(bbox['right'], bbox['bottom']),
(bbox['right'], bbox['top']),
(bbox['left'], bbox['top']),
(bbox['left'], bbox['bottom']),
))
extended_grid_extent = Polygon(((extended_bbox['left'], extended_bbox['bottom']),
(extended_bbox['right'], extended_bbox['bottom']),
(extended_bbox['right'], extended_bbox['top']),
(extended_bbox['left'], extended_bbox['top']),
(extended_bbox['left'], extended_bbox['bottom']),
))
In [13]:
rows
Out[13]:
In [14]:
columns
Out[14]:
In [15]:
extended_rows
Out[15]:
In [16]:
extended_columns
Out[16]:
In [ ]:
#values = pandas.DataFrame({'id':numpy.arange(rows*columns), 'naturalness_500m':numpy.zeros(rows*columns)}).set_index('id')
values = pandas.DataFrame({'id':numpy.arange(extended_rows * extended_columns), 'naturalness_500m':numpy.zeros(extended_rows * extended_columns)}).set_index('id')
In [ ]:
n_calculations = extended_rows * extended_columns
progress = pyprind.ProgBar(n_calculations, bar_char='█', title='progress', monitor=True, stream=1, width=70)
for i, coords in enumerate(extended_grid.ravel()):
progress.update(item_id=i)
point = Point(coords)
#if point.within(grid_extent):
# value = 999
#else:
value = get_naturalness(point=point, radius=radius)
values['naturalness_500m'].iloc[i] = value
In [ ]:
# save values
values.to_csv("/home/ubuntu/data/model grid/naturalness_500m_grid_200m_spacing.csv")
In [12]:
data = values['naturalness_500m'].as_matrix().reshape((rows, columns))
data = numpy.flipud(data)
In [13]:
plt1 = pyplot.imshow(data, cmap='viridis')
clb1 = pyplot.colorbar()
In [19]:
# import values
values = pandas.read_csv("/home/ubuntu/data/model grid/naturalness_500m_grid_200m_spacing_merged.csv")
In [20]:
# reshape into grid
data = values['naturalness_500m'].as_matrix().reshape((extended_rows, extended_columns))
# mean center values (mean computed in the regression model notebook)
u = data - 2.647528316787191
# define model parameters (computed in the regresssion model notebook)
g_00_b = -0.4061 # level 2
g_01_b = 1.2439
g_10_b = 0.1696
g_11_b = 0.2671
a_b = g_00_b + (g_01_b * u)
b_b = g_10_b + (g_11_b * u)
x = [t for t in numpy.linspace(-3, 3, 21)] # level 1
y_b = numpy.empty(shape=(len(x), u.shape[0], u.shape[1]), dtype=numpy.float64)
for i, xi in enumerate(x):
y_b[i] = a_b + (b_b * xi)
export_to_hdf5(filepath="/home/ubuntu/predicted_biophony.hd5", name="predicted biophony", data=y_b)
transform predicted biophony to range from 0 to 1
In [21]:
y_b = y_b - y_b.min()
y_b = y_b / y_b.max()
define folium map
In [22]:
# center the map on the site named 'Hofgarten' and use the "Stamen Toner" background tiles
map_center = Site.objects.get(name='Hofgarten').geometry.coords
# transfrom the points from MGI / Austria GK West to WGS 84
lon, lat = pyproj.transform(p1, p2, map_center[0], map_center[1])
# intialize the map and use the "Stamen Toner" background tiles
map_b = folium.Map(location=[lat, lon],
zoom_start=15, tiles="Stamen Toner",
detect_retina=True)
create the biophony overlay (for x[0]) and add it to the map
In [23]:
b_overlay = ImageOverlay(y_b[15],
bounds=[numpy.roll(numpy.array(\
pyproj.transform(p1, p2, extended_bbox['left'], extended_bbox['bottom'])), 1).tolist(),
numpy.roll(numpy.array(\
pyproj.transform(p1, p2, extended_bbox['right'], extended_bbox['top'])), 1).tolist()],
opacity=0.5,
origin='lower',
colormap=cm.Greens).add_to(map_b)
show the map
In [24]:
map_b
Out[24]:
create and save images in PNG format for each time step (each value of x)
In [ ]:
image_path = ""
for i in range(y.shape[0]):
image = Image.fromarray(cm.Greens(numpy.flipud(y[i]), bytes=True))
image.save(path.join(image_path, "image_{0}.png".format(i + 1)))
In [25]:
# import values
values = pandas.read_csv("/home/ubuntu/data/model grid/pavement_100m_grid_200m_spacing_merged.csv")
In [26]:
# reshape into grid
data = values['pavement_100m'].as_matrix().reshape((extended_rows, extended_columns))
# mean center values (mean computed in the regression model notebook)
u = data - (-48.8641)
# define model parameters (computed in the regresssion model notebook)
g_00_a = -0.0876 # level 2
g_01_a = -0.3314
y_a = g_00_a + (g_01_a * u) # level 1
export_to_hdf5(filepath="/home/ubuntu/sel.hd5", name="predicted SEL", data=y_a)
transform predicted biophony to range from 0 to 1
In [27]:
# transform to range 0-1
y_a = y_a + abs(y_a.min())
y_a = y_a / y_a.max()
# transform to account for scale
#y_a = (y_a * 0.50) + 0.3
define folium map
In [28]:
# center the map on the site named 'Hofgarten' and use the "Stamen Toner" background tiles
map_center = Site.objects.get(name='Hofgarten').geometry.coords
# transfrom the points from MGI / Austria GK West to WGS 84
lon, lat = pyproj.transform(p1, p2, map_center[0], map_center[1])
# intialize the map and use the "Stamen Toner" background tiles
map_a = folium.Map(location=[lat, lon],
zoom_start=15, tiles="Stamen Toner",
detect_retina=True)
create the anthrophony overlay and add it to the map
In [29]:
a_overlay = ImageOverlay(y_a,
bounds=[numpy.roll(numpy.array(\
pyproj.transform(p1, p2, extended_bbox['left'], extended_bbox['bottom'])), 1).tolist(),
numpy.roll(numpy.array(\
pyproj.transform(p1, p2, extended_bbox['right'], extended_bbox['top'])), 1).tolist()],
opacity=0.5,
origin='lower',
colormap=cm.YlOrBr).add_to(map_a)
show the map
In [30]:
map_a
Out[30]:
create and save an image in PNG format of the overlay
In [ ]:
image = Image.fromarray(cm.YlOrBr(numpy.flipud(y_a), bytes=True))
image.save(path.join(image_path, "image_1_a.png"))