# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: 'Python 3.6.13 64-bit (''fossworkshop'': conda)'
#     name: python3
# ---

# # Vector Data Analysis
# ## This notebook represents working with vector data in python
#
# Vector data is usually a tabular data coupled with location information. e.g. Data of all states in India ( This file will have some attribute data about states such as name, population, etc. along with one column of geometry containing location information). Vector data geometry can be divided in 3 major types: 
#
# 1. Point Geometry  
#
#     Point Geometry consists of discrete location information such as *latitude, longitude* which can help us to identify the exact location of given feature. 
#
#     e.g - Location of bus stop, location of user, etc.
#
# 2. Line Geometry
#
#     Line Geometry is collection of multiple *latitude, longitude* in an array which represents continuous path.
#
#     e.g. - Centreline of road, River, Path created by user, etc.
#
# 3. Polygon Geometry 
#
#     Polygon Geometry is collection of multiple *latitude, longitude* in an array which represents continuous enclosed area.
#
#     e.g. - Geometry of State, polygon of building, etc.
#     
# ## Loading Data
#
# First step is to load the Data in to python, this data can be a file available on machine, data stored in database, or file hosted on some server
#
# ### 1. Loading Shapefile 
#
# Loading all countries geometry (src: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/)

#Using geopandas to work with data
import geopandas as gpd
#Using Matplotlib for visualisation
import matplotlib
# %matplotlib inline


# +
#load it as a pandas dataframe with understanding on geometrical data
countries = gpd.read_file('../data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp')

countries
# -

countries.plot()


countries.head()

countries.tail()

# #### Understanding GeoDataFrame
#
# GeoDataFrame will always have <b>geometry</b> column, apart from that other columns will act as metadata. 
#
# So `geopandas` = `pandas` + `geometry` 
#
# Each column except geometry in the geopandas is of type `pandas.Series` , geometry is treated as `pandas.GeoSeries`

print(type(countries.geometry))
print(type(countries.scalerank))


# Each geometry is a `shapely` Shape, thus we can perform all shapely methods on these geometries
#
# Checkout all available methods here https://shapely.readthedocs.io/en/stable/manual.html#predicates-and-relationships

countries.geometry

countries.geometry.centroid

# ### 2. Loading Geojson 
#
# Loading local geojson file 

rivers = gpd.read_file('../data/rivers.geojson')
rivers

# ### 3. Loading PostgreSQL
#
# Loading data from database

# +
import psycopg2 

con = psycopg2.connect(database="postgres", user="postgres", password="<PASSWORD>",
    host="localhost")

sql = "SELECT * FROM public.places"
places = gpd.read_postgis(sql, con, geom_col='geom' )
places
# -

# ### 4. Import CSV 
#
# Assuming that CSV has a geometry column that contains geometery in WKT format

# +
from shapely import wkt

airport = gpd.read_file('../data/airport.csv')

airport['geometry'] = airport['geom'].apply(wkt.loads)
del airport['geom']
airport

# -

# ### 5. Creating geometry On the fly
#
# Create geodataframe from csv having columns as longitude and latitude, which will be used further to create geometery on the fly

import pandas as pd
df = pd.read_csv('../data/stadium.csv')
stadium = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df.lon, df.lat))
stadium

# ### 6. Create Geodataframe manually 
#
# User can also create Geodataframe in the notebook, using their own data

# +
from shapely.geometry import Point

police = gpd.GeoDataFrame({
    'geometry': [Point(1, 1), Point(2, 2),Point(2, 1),Point(1, 2),Point(1.5, 2)],
    'id': [1, 2,3,4,5],
    'criminals': [12,34,112,41, 212]})
police
# -

police.to_html('police.html')

# ## More about shapely
#
# ### How to create geometery



# +
from shapely.geometry import Polygon,Point,LineString
Pt = Point(10,10)
line = LineString([(0,0),(0,3),(3,0)])
poly = Polygon([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])


# -

Pt

line

poly

# ### Geospatial analysis

poly.touches(line)

poly.contains(Pt)

Pt.buffer(20).contains(poly)

# ### Peoperties of shape
#

poly.area

line.bounds

line.length

# ## More about Fiona 
#
# Fiona is a python interface of GDAL/OGR library, Geopandas is a more easy to user wrapper

import fiona

places = fiona.open('../data/ne_10m_populated_places/ne_10m_populated_places.shp')
places

places.driver

places.schema

places.crs

len(places)

# ## Playing with GeoDataFrame
#
# ### Coordinate system 
#
# Unline `shapely`, `geopandas` understands crs
#
# What CRS are important? 
# - CRS will make sense out of your data such as whether the units are degrees/meters
# - Bringing all data in same CRS allows us to do spatial analysis with data 
#
#
# To check CRS of GeoDataframe
#

countries.crs

# We can also set CRS for the GeoDataFrame which has no default CRS

We can also set CRS for the GeoDataFrame which has no default CRS

police = police.set_crs('epsg:4326')
police.crs

# We can also convert GeoDataFrame from one CRS to another

police.crs

police_3857 = police.to_crs(3857)

police_3857

police_3857.crs

# ### Merging 
#
# 1. Atrribute based merge

neighbor = pd.DataFrame({
    'id': [1, 2,3,4,5],
    'neighbor_id': ['a1', 'a2','a3','c4','d5'],
    'neighbor_name': ['andy','julio','true','skewd', 'tauras']})
neighbor

updated_police = police.merge(neighbor, on='id')
updated_police

# 2. Spatial merge

pd.set_option('max_columns', 100)
airport = airport.set_crs('epsg:4326')
airport.head()

simple_countries = countries[['ADMIN','geometry']]
simple_countries.head()

# +
airport_with_country = gpd.sjoin(airport, simple_countries, how="inner", op='intersects')


# -

airport_with_country.head()

# <b>op</b> : Another way to perform same query can be using operation `within` instead of `intersect` .

airport_with_country_within = gpd.sjoin(airport, simple_countries, how="inner", op='within')


airport_with_country_within.tail()

# <b>how</b> : We can use `left` , `right` , `inner` .
#
#  `left`: use the index from the first (or left_df) geodataframe that you provide to sjoin; retain only the left_df geometry column
#
# `right`: use index from second (or right_df); retain only the right_df geometry column
#
# `inner`: use intersection of index values from both geodataframes; retain only the left_df geometry column

airport_with_country_right = gpd.sjoin(airport, simple_countries, how="right", op='within')
airport_with_country_right.head()


# ### Edit the existing data
#
# Editing metadata

updated_police.iloc[0]

updated_police.iloc[0,2] = 24

updated_police.iloc[0]

# Editing geometry

# +
from shapely.geometry import Point

updated_point = Point(3,4)
updated_police.iloc[0,0] = updated_point
updated_police
# -

# ### Querying data
#
# 1. Based on metadata

countries.head()

India = countries[countries['ADMIN'] == "India"]
India

densly_pop = countries[countries['POP_EST'] > 100000000]
densly_pop

countriesWithC = countries[countries['SOVEREIGNT'].str.startswith('C')]
countriesWithC

densecountriesWithC = countries[(countries['SOVEREIGNT'].str.startswith('C')) &  (countries['POP_EST'] > 1000000000)]
densecountriesWithC

# 2. Spatial Query
#
# Spatial query uses shapely geometry as base geometry on top of which geodataframe can be queried.
# Available oprations are listed at
# https://shapely.readthedocs.io/en/latest/manual.html#binary-predicates

indian_shape = India['geometry'].squeeze()

type(India['geometry'].squeeze())

test_pt = Point(1,1)

test_pt.intersects(indian_shape)

nashik = Point(73.76,19.96)

nashik.within(indian_shape)

indian_airport = airport[airport.within(indian_shape)]

indian_airport

# #### Quiz -> Can you create the dataframe of all airports and cities within your country

indian_river = rivers[rivers.intersects(indian_shape)]
indian_river.plot()

Neighbours_India = countries[countries.touches(indian_shape)]
Neighbours_India.plot()

# ### Geospatial Operations
#
# Understanding base logic first! Back to `shapely`

test_point = Point(0,0)
test_point

test_point.buffer(10)

test_point.buffer(10).area

# +
from shapely.geometry import LineString

test_line = LineString([(0, 0), (1, 1), (0, 2)])
test_line
# -

#Buffer puts original geometry at center and create buffer alongside
test_line.buffer(0.1)

# +
#We can also put geometry on either side ( Positive value will put buffer to left)

test_line.buffer(0.5, single_sided=True)

# +
#We can also put geometry on either side ( negative value will put buffer to right)

test_line.buffer(-0.5, single_sided=True)
# -

# Operations on `geopandas`

Indian_cities =  places[places.within(indian_shape)]

Indian_cities

Indian_cities_m = Indian_cities.to_crs(3857)

Indian_cities_m.crs

Indian_cities_m.head()

city_buffer = Indian_cities_m[['geom','name']]
city_buffer

fig, ax = plt.subplots(figsize=(16, 16))
India_m.plot(ax=ax, color='#ffffff', edgecolor='#6a6a6a', linewidth=2)
city_buffer.plot(ax=ax, color='#f00', edgecolor='#000000')


city_buffer['geom'] = city_buffer.buffer(50000)

city_buffer

countries.head()

countries_centroid = countries[['geometry','NAME','CONTINENT']]
countries_centroid.head()

countries_centroid['geometry'] = countries_centroid['geometry'].centroid
countries_centroid.head()

countries_centroid.plot()

countries['area'] = countries['geometry'].area
countries.head()

countries_m = countries.to_crs(3857)
countries_m['area'] = (countries_m['geometry'].area)/1000000
countries_m

# ## Visualising GeoDataFrame

#simple visualisation 
countries_m.plot()

countries_m = countries_m[countries_m['NAME'] != "Antarctica"]
countries_m.plot()

#color based on column
countries_m.plot(column='CONTINENT')

countries_m.plot(column='CONTINENT',legend=True)

# +
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(16, 16))
countries_m.plot(ax=ax,column='CONTINENT',legend=True)
# -

ax = countries_m.plot(column='CONTINENT',legend=True)
ax.set_axis_off()

#Checkout available color maps => https://matplotlib.org/2.0.2/users/colormaps.html
countries_m.plot(column='CONTINENT',  cmap='winter')


countries_m.plot(column='POP_EST',legend=True)


countries_plot = countries_m[(countries_m['NAME'] != 'India') & (countries_m['NAME'] != 'China')]

countries_plot.plot(column='POP_EST',legend=True,figsize=(16,16), legend_kwds={'label': 'Population'})

# ### matplotlib to show multiple data 

basemap = countries_m.plot(column='CONTINENT', cmap='cool')
cities_m = places.to_crs(3857)
cities_m.plot(ax=basemap, marker='o', color='red', markersize=5)

#load world polygon
bbox = gpd.read_file('../data/world.geojson')
world = bbox.loc[0].geometry
world

cities_m = cities_m[cities_m.within(world)]

basemap = countries_m.plot(column='CONTINENT', cmap='cool')
cities_m.plot(ax=basemap, marker='o', color='red', markersize=5)

# ### geopandas overlay to show multiple data 

# +
fig, ax = plt.subplots(figsize=(16, 16))

India_m.plot(ax=ax, color='b', edgecolor='#f0f', linewidth=2)
Indian_cities_m.plot(ax=ax, color='r', edgecolor='#fff')

# -

Indian_cities_m['geom'] = Indian_cities_m['geom'].buffer(50000)

# +
fig, ax = plt.subplots(figsize=(16, 16))

India_m.plot(ax=ax, color='b', edgecolor='#f0f', linewidth=2)
Indian_cities_m.plot(ax=ax, color='r', edgecolor='#fff')
# -

non_rural_area = gpd.overlay(India_m, Indian_cities_m, how='difference')
non_rural_area.plot(figsize=(16, 16))

# ## Interactive Maps in python

from ipyleaflet import Map, GeoData, basemaps, LayersControl
import geopandas as gpd 


# ## Loading empty Map
#
# This will include initiating map with center, zoom level and basemap choice. 
# Checkout available basemap options at : https://ipyleaflet.readthedocs.io/en/latest/api_reference/basemaps.html?highlight=basemap

m = Map(center=(27,71), zoom = 3, basemap= basemaps.Esri.WorldTopoMap)
m

# ## Loading Data to Map
#
# ### 1. Loading Geopandas dataframe

# +

m = Map(center=(27,71), zoom = 3, basemap= basemaps.Stamen.Toner)
cities = gpd.read_file('../data/ne_10m_populated_places/ne_10m_populated_places.shp')
cities_data  = GeoData(geo_dataframe = cities,
    style={'color': 'black', 'radius':4, 'fillColor': '#3366cc', 'opacity':0.5, 'weight':1.9, 'dashArray':'2', 'fillOpacity':0.6},
    hover_style={'fillColor': 'red' , 'fillOpacity': 0.2},
    point_style={'radius': 5, 'color': 'red', 'fillOpacity': 0.8, 'fillColor': 'blue', 'weight': 3},
    name = 'Release')
m.add_layer(cities_data)
m
# -

# ### 2. Loading WMS layer
#

# +
from ipyleaflet import Map, WMSLayer, basemaps

wms = WMSLayer(
    url='https://ahocevar.com/geoserver/wms',
    layers='topp:states',
    format='image/png',
    transparent=True,
    attribution='Made for GeoPython 2021'
)

m = Map(basemap=basemaps.CartoDB.Positron, center=(38.491, -95.712), zoom=4)

m.add_layer(wms)

m
# -

# ## Adding Popup 
#
# ### 1. Adding static popup

from ipywidgets import HTML
from ipyleaflet import Map, Marker, Popup
center = (19.975040, 73.763190)
m = Map(center=center, zoom=17, close_popup_on_click=False)
marker = Marker(location=(19.975040, 73.763190))
m.add_layer(marker)
message2 = HTML()
message2.value = "Hey!! I'm speaking at foss4g 2021 🔥"
marker.popup = message2
m

# ### 2. Using Custom data for popup
#
# For this example we'll prepare map of following scenario
# Seeing all the cities as a point on map and on click show their name

# +
#Preparing data 
all_cities = gpd.read_file('../data/ne_10m_populated_places/ne_10m_populated_places.shp')
all_countries =  gpd.read_file('../data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp')
all_cities.dropna(subset=["NAME","geometry"])
India = all_countries[all_countries['NAME'] == 'India']
Indian_cities = all_cities[all_cities.within(India.squeeze().geometry)]
Indian_cities

# Creating Map
from ipyleaflet import Map, Marker, Popup
from ipywidgets import HTML
center = (33.762918,68.637469)
m = Map(center=center, zoom=3, close_popup_on_click=False)

# Adding data as marker 
for index, row in Indian_cities.iterrows():
    message2 = HTML()
    marker = Marker(location=(row['geometry'].y, row['geometry'].x))
    message2.value = row['NAME']
    # message2.description = row['NAME']
    marker.popup = message2
    m.add_layer(marker)
#     print(index)

#load map
m
# -

# ## Another interesting map options
#
# 1. AntPath 
# 2. Marker Cluster
# 3. Heatmap
# 4. Velocity
# 5. Choropleth
#
# check out out at https://ipyleaflet.readthedocs.io/
#
# ## Controls in map
#
# Different Controls can be added to the map to make it more user friendly. Some of such controls are as follows
#
# ### 1. Scale control

# +
from ipyleaflet import Map, ScaleControl

m = Map(zoom=15, center=[19.975040, 73.763190])
m.add_control(ScaleControl(position='bottomleft'))

m
# -

# ### 2. Split Map

# +
from ipyleaflet import Map, basemaps, basemap_to_tiles, SplitMapControl

m = Map(zoom=15, center=[19.975040, 73.763190])

right_layer = basemap_to_tiles( basemaps.Stamen.Toner)
left_layer = basemap_to_tiles(basemaps.CartoDB.Positron)

control = SplitMapControl(left_layer=left_layer, right_layer=right_layer)
m.add_control(control)

m
# -

# ### Apart from these, some of the most widely used controls are
#
# 1. Draw on map
# 2. Adding Legends
# 3. Measure, etc. 
#
# You can find all available controls at https://ipyleaflet.readthedocs.io/en/latest/index.html (Look for control section)

# ## pydeck
#
# Python package based on deck.gl https://pydeck.gl/ which also provides support for 3d data and visualisation
#
# install the package `pip install pydeck`
#
# pydeck by default uses carto basemap, but it can be replaced with `Mapbox` or `Google`, to do so, you will need to get API key from their website

# +

import pydeck as pdk
import pandas as pd
data = '../data/flights.csv'
commute_pattern = pd.read_csv(data)


# view (location, zoom level, etc.)
view = pdk.ViewState(latitude=21.214885, longitude=77.950061, pitch=50, zoom=3)

# layer
# from home (orange) to work (purple)
arc_layer = pdk.Layer('ArcLayer',
                      data=commute_pattern,
                      get_source_position=['lon_from', 'lat_from'],
                      get_target_position=['lon_to', 'lat_to'],
                      get_width=5,
                      get_tilt=15,
                      pickable=True,
                         auto_highlight=True,
                      # RGBA colors (red, green, blue, alpha)
                      get_source_color=[255, 165, 0, 80],
                      get_target_color=[128, 0, 128, 80])

# render map
# choose map style
TOOLTIP_TEXT = {"html": "{flights} flights taken <br /> on this root"}

arc_layer_map = pdk.Deck(
                         layers=arc_layer,
                         initial_view_state=view,
                         tooltip=TOOLTIP_TEXT
                    )
arc_layer_map.to_html('deck.html')
arc_layer_map.show()
# -

# checkout other options at : https://pydeck.gl/index.html



StarcoderdataJupyterScriptsDedupFiltered