Create a gist now

Instantly share code, notes, and snippets.

Embed
Grab information on forest cover loss for the supplied shapefile
import os
import requests
import json
import pandas as pd
# Go from JSON to a measurement of UMD deforestation by year. If you
# have a shapefile, you will have to convert and simplify the
# polygons, preserving topology. The Earth Engine API does not
# support long, long requests. You can convert the shapefile using
# OGR within the data subdirectory:
# ogr2ogr -f GeoJSON map.geojson map.shp -simplify 0.001
# You can also upload your various formats (csv, shp, kml, etc.) to
# CartoDB and export to JSON. Regardless of the path you choose,
# young Padawan, JSON is the destination. Specifically,
# data/map.json.
def _read_data(data='data/map.geojson'):
"""Returns a list of JSON dictionaries, one for each feature of
the polygons.
"""
with open(data) as json_file:
x = json.load(json_file)
return x['features']
def _force_ring(coords):
"""Returns a nested list of coordinates, ensuring that the
coordinates form a ring with the first and last coordinates the
same.
"""
if coords[0] != coords[-1]:
return [coords + coords[-1]]
else:
return [coords]
def _polygon(coords):
"""Converts a list of coordinates to json"""
return json.dumps({"type": "Polygon", "coordinates": _force_ring(coords)})
def _params(coords, year):
"""Returns a parameter dictionary for a post request"""
p = _polygon(coords)
return {"begin":year-1, "end":year, "geom":p}
def _grab_loss(coords, year):
"""Returns the loss in hectares of the supplied coordinates and
year from the UMD data set hosted on Earth Engine via GFW API.
"""
endpoint = 'http://gfw-apis.appspot.com/datasets/umd'
res = requests.post(endpoint, data=_params(coords, year))
return res.json()['loss']
def _process_entry(entry):
"""Converts the supplied entry, a multipolygon in json, into a
list of dictionaries with loss information associated with year
and province identifying information. Ready for quick conversion
to Pandas data frame.
"""
name_1 = entry['properties']['NAME_1']
name_2 = entry['properties']['NAME_2']
coord_set = entry['geometry']['coordinates']
res = []
for yr in range(2001,2013):
x = [_grab_loss(coord, yr) for coord in coord_set]
xx = {'prov':name_1, 'subprov':name_2, 'year':yr, 'loss':sum(x)}
res.append(xx)
return res
def _filter_admin(prov, subprov, data='data/map.geojson'):
"""Accepts a province and subprovince name (strings) and a data
source in json format and returns the multipolygon associated with
that administrative unit.
"""
polys = _read_data(data)
def _spec_filter(xx):
x = xx['properties']
return (x['NAME_1'] == prov) & (x['NAME_2'] == subprov)
return filter(_spec_filter, polys)[0]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment