import requests
import json
import pandas as pd
import geopandas as gpd
import shapely
from shapely.geometry import shape, Point
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import holoviews as hv
import iris
import cartopy
from cartopy import crs
from cartopy import feature as cf
from bokeh.models import WMTSTileSource
from bokeh.embed import file_html
from bokeh.plotting import output_file
from bokeh.io import save
import shapefile
from IPython.display import HTML
import warnings
warnings.filterwarnings('ignore')
# This function fetches a shapefile from the Bushtel site, and processes it to make it a GeoPandas dataframe
def get_community_geodataframe(community_id):
url = "http://www.bushtel.nt.gov.au/api/Lot?communityIds={0}".format(str(community_id))
r = requests.get(url)
j = json.loads(r.content.decode())
for item in j:
item['Bounds'] = json.loads(item['Bounds'])
for item in j:
item['geometry'] = shape(item['Bounds'])
for item in j:
item['geom_type'] = str(item['geometry'])
gdf = gpd.GeoDataFrame(j)
gdf.LotNumber = gdf.LotNumber.astype(int)
gdf.CommunityTypeId = gdf.CommunityTypeId.astype(int)
gdf = gdf[[ 'Bounds', 'CommunityName', 'CommunityTypeId', 'CommunityTypeName', 'Id', 'LaisCode', 'LotNumber', 'LotTypeId', 'LtoCode', 'StreetAddress', 'geom_type', 'geometry']]
gdf = gdf[(gdf.geom_type == 'Polygon')]
#gdf.crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
gdf.crs = '+init=epsg:4326'
return gdf
# Get the map for Gapuwiyak
bt_gdf = get_community_geodataframe(500)
# Trim out road polygon
bt_gdf = bt_gdf[(bt_gdf.index != 96)]
# Fetch an exported CSV of the data from the GeoODK App
df = pd.read_csv('GAPUWIYAK_results.csv', header=0)
odk_gdf = gpd.GeoDataFrame(df)
#print(list(odk_gdf))
odk_gdf['geometry'] = [Point(xy) for xy in zip(odk_gdf['coordinates:Longitude'], odk_gdf['coordinates:Latitude'])]
odk_gdf.crs = '+init=epsg:4326'
# Merge the data from the Bushtel site with the data from the GeoODK App
gdf = gpd.sjoin(bt_gdf, odk_gdf, how="left", op='intersects')
# Create Nice Names for Fields
gdf['Lot Number'] = gdf[['LotNumber']]
gdf['Building Colour'] = gdf[['build_colour']]
gdf['Points of Interest'] = gdf[['points_interest_other']]
def changeColor(df):
if (df[['build_colour']] == "d_blue").all():
return "Dark Blue"
if (df[['build_colour']] == "red").all():
return "Red"
if (df[['build_colour']] == "yellow").all():
return "Yellow"
if (df[['build_colour']] == "l_green").all():
return "Light Green"
if (df[['build_colour']] == "orange").all():
return "Orange"
if (df[['build_colour']] == "d_green").all():
return "Dark Green"
if (df[['build_colour']] == "grey").all():
return "Grey"
if (df[['build_colour']] == "l_blue").all():
return "Light Blue"
else:
return "None"
gdf['Colour'] = gdf.apply(changeColor,axis=1)
gdf = gdf.fillna('None')
# Create shapefiles for each of the colours of house
gdf_red = gdf[(gdf.build_colour == 'red')]
gdf_l_blue = gdf[(gdf.build_colour == 'l_blue')]
gdf_d_blue = gdf[(gdf.build_colour == 'd_blue')]
gdf_l_green = gdf[(gdf.build_colour == 'l_green')]
gdf_d_green = gdf[(gdf.build_colour == 'd_green')]
gdf_yellow = gdf[(gdf.build_colour == 'yellow')]
gdf_brown = gdf[(gdf.build_colour == 'brown')]
gdf_grey = gdf[(gdf.build_colour == 'grey')]
gdf_white = gdf[(gdf.build_colour == 'white')]
gdf_orange = gdf[(gdf.build_colour == 'orange')]
gdf_other = gdf[(gdf.build_colour != 'red') &
(gdf.build_colour != 'l_blue') &
(gdf.build_colour != 'd_blue') &
(gdf.build_colour != 'l_green') &
(gdf.build_colour != 'd_green') &
(gdf.build_colour != 'yellow') &
(gdf.build_colour != 'brown') &
(gdf.build_colour != 'grey') &
(gdf.build_colour != 'white') &
(gdf.build_colour != 'orange')]
# Export the shapefiles
gdf_red.to_file('red_shape')
gdf_l_blue.to_file('l_blue_shape')
gdf_d_blue.to_file('d_blue_shape')
gdf_l_green.to_file('l_green_shape')
gdf_d_green.to_file('d_green_shape')
gdf_yellow.to_file('yellow_shape')
#gdf_brown.to_file('brown_shape')
gdf_grey.to_file('grey_shape')
#gdf_white.to_file('white_shape')
gdf_orange.to_file('orange_shape')
gdf_other.to_file('other_shape')
# Set up the Esri Worldmap tile source
ts = WMTSTileSource(url='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{Z}/{Y}/{X}.jpg')
# Set up GeoViews
import geoviews as gv
hv.notebook_extension('bokeh', logo=False)
%output backend='bokeh'
# Create main overlay for tooltips
sf = cartopy.io.shapereader.Reader('join_gap_shape/join_gap_shape.shp')
tooltip_sfpolys = gv.Shape.from_records(sf.records(), gdf, on='Id', value='LotNumber', index=['Lot Number', 'Colour', 'Points of Interest'], crs=crs.PlateCarree())
tooltip_sfpolys = tooltip_sfpolys.opts({'Shape': {'style': {'fill_color':'red', 'fill_alpha':0.0}, 'plot':dict(tools=['hover'])} })
# Function to create polygon layers for each colour
def create_sfpoly(shapefile, color, alpha=0.4):
reader = cartopy.io.shapereader.Reader(shapefile + '/' + shapefile + '.shp')
poly_diag = gv.Shape.from_records(reader.records(), crs=crs.PlateCarree())
poly_diag = poly_diag.opts({'Shape': {'style': {'fill_color':color, 'fill_alpha':alpha}} })
return poly_diag
# Create all the polygon layers for the houses
red_sfpolys = create_sfpoly(shapefile='red_shape', color='red')
l_blue_sfpolys = create_sfpoly(shapefile='l_blue_shape', color='cyan')
d_blue_sfpolys = create_sfpoly(shapefile='d_blue_shape', color='blue')
d_green_sfpolys = create_sfpoly(shapefile='d_green_shape', color='darkgreen')
l_green_sfpolys = create_sfpoly(shapefile='l_green_shape', color='lawngreen')
grey_sfpolys = create_sfpoly(shapefile='grey_shape', color='grey', alpha=0.6)
orange_sfpolys = create_sfpoly(shapefile='orange_shape', color='orange')
yellow_sfpolys = create_sfpoly(shapefile='yellow_shape', color='yellow')