This notebook compares the land cover within the entire study area to the land cover within a certain radius around sample locations.
In [24]:
from database.models import Site
from landscape.models import LandCover
from geo.models import LandCoverType
from geo.models import Boundary
from django.db import connection
from django.contrib.gis.geos import GEOSGeometry
from os import path
import numpy
import pandas
import geopandas
import fiona
import pyproj
from matplotlib import pyplot
%matplotlib inline
define the project coordinate system
In [2]:
austria_mgd = "+init=epsg:31254"
define the output filepath
In [3]:
output_filepath = "/Users/Jake/OneDrive/Documents/alpine soundscapes/data/landcover"
define schema for GeoJSON output
In [15]:
schema = {'geometry': "MultiPolygon",
'properties': {'id': 'int',
'cover_type': 'int'}}
In [5]:
# returns a geopandas dataframe from a given PostGIS query
def get_geodataframe(queryset, modification=None, index_col=None, crs=austria_mgd):
query = queryset.query.sql_with_params()
if modification:
query = (modification, query[1])
return geopandas.read_postgis(query[0], connection,
geom_col='geometry',
params=query[1],
index_col=index_col,
crs=crs)
In [6]:
# returns a geopandas dataframe containing the land cover intersected with a given geometry
def get_landcover(geometry, make_valid=False, index_col=None, crs=austria_mgd):
if make_valid:
intersection = """ST_Intersection(ST_MakeValid("landscape_landcover"."geometry"), %s)"""
else:
intersection = """ST_Intersection("landscape_landcover"."geometry", %s)"""
query = """
SELECT "landscape_landcover"."id", "landscape_landcover"."cover_type",
{0} AS "geometry" FROM "landscape_landcover"
WHERE ST_Intersects("landscape_landcover"."geometry", %s)
""".replace('\n', '').format(intersection)
return geopandas.read_postgis(query, connection,
geom_col='geometry',
params=["SRID={0};{1}".format(geometry.srid, geometry.wkt),
"SRID={0};{1}".format(geometry.srid, geometry.wkt)],
index_col=index_col,
crs=crs)
In [31]:
# returns a numpy array with the percentage of each land cover type in respect to the total geometry area
def get_percentages(source):
y = range(1, 16)
x = numpy.empty(15)
total_area = source.area.sum()
for i, y_i in enumerate(y):
x[i] = (source[source.cover_type == y_i].area.sum() / total_area) * 100
return x
define a list of radii defining the surrounding land cover
In [7]:
radii = [50, 100, 200, 500]
compute a geometry (unioned buffer around locations) for each radius
In [8]:
location_area_list = []
for radius in radii:
location_areas = get_geodataframe(Site.objects.filter(id__lte=30)).buffer(radius)
location_area = location_areas.iloc[0]
for area in location_areas:
location_area = location_area.union(area)
# convert to geos multipolygon with srid
location_area = GEOSGeometry(location_area.wkt)
location_area.srid = 31254
location_area_list.append(location_area)
query database
In [34]:
location_landcover_list = []
for location_area, radius in zip(location_area_list, radii):
location_landcover = get_landcover(location_area)
location_landcover_list.append(location_landcover)
# export result to a GeoJSON file
#location_landcover.to_file(filename=path.join(output_filepath, "location_landcover_{0}m.geojson".format(radius)),
# driver='GeoJSON', schema=schema)
get the study area boundary geometry from the database
In [10]:
boundary = Boundary.objects.get(name="study area")
query database
In [16]:
study_landcover = get_landcover(boundary.geometry, make_valid=True)
# export result to a GeoJSON file
#study_landcover.to_file(filename=path.join(output_filepath, "study_area_landcover.geojson".format(radius)),
# driver='GeoJSON', schema=schema)
In [32]:
area_percentages = pandas.DataFrame({'name': [r['name'].lower() for r in LandCoverType.objects.filter(id__lte=15).values()],
'study_area': get_percentages(study_landcover)}, index=range(1, 16))
In [35]:
for i, radius in enumerate(radii):
area_percentages["{0}m".format(radius)] = get_percentages(location_landcover_list[i])
In [36]:
area_percentages
Out[36]:
In [37]:
area_percentages.to_csv(path.join(output_filepath, "area_percentages.csv"))