Map the acoustic environment

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.

Import statements


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

Variable declarations

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

Function declarations


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()

Create grid of points

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)]

extended grid points


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]:
47

In [14]:
columns


Out[14]:
58

In [15]:
extended_rows


Out[15]:
60

In [16]:
extended_columns


Out[16]:
75

Compute land cover area for all points in the grid


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


progress
0%                                                                  100%
[███                                                                   ] | ETA: 00:03:53 | Item ID: 192

In [ ]:
# save values
values.to_csv("/home/ubuntu/data/model grid/naturalness_500m_grid_200m_spacing.csv")

View result


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()


Predict biophony


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)

Map predicted biophony

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)))

Predict anthrophony


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)

Map anthrophony

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"))