import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

%matplotlib inline

districts = gpd.read_file("Community Districts/districts.shp")

districts.head()

districts[districts.boro_cd > 400]

districts.head()

districts.plot()

# Go into the crs to convert it...
# ignore the datum and spheroid,
# just change the PROJECTION to MERCATOR
districts.to_crs({'proj': 'merc'}).plot()

# Give it the SECRET CODE from the PETROLEUM GROUP
# (which you can try to find by googling)
# (or hopefully you have a list because they're
# all very confusingly/similarly named)
districts.to_crs(epsg=3857).plot()

df = pd.read_csv("../earthquakes_1.0_day.csv")
df.head(2)

def make_point(row):
    return Point(row.longitude, row.latitude)

# Go through every row, and make a point out of its lat and lon
points = df.apply(make_point, axis=1)

# Make a new GeoDataFrame
# using the data from our old df
# but also adding in the geometry we just made
earthquakes = gpd.GeoDataFrame(df, geometry=points)

# It doesn't come with a CRS because it's a CSV, so let's
# say "hey, let's use the standard shape of the earth etc"
earthquakes.crs = {'init': 'epsg:4326'}

# Let's look at the first few
earthquakes.head()

earthquakes.plot()

# Read in the CSV
df = pd.read_csv("../earthquakes_1.0_day.csv")

points = df.apply(lambda row: Point(row.longitude, row.latitude), axis=1)
earthquakes = gpd.GeoDataFrame(df, geometry=points)
earthquakes.crs = {'init': 'epsg:4326'}

# If you want to know how this all works, look above
earthquakes.head(2)

earthquakes.plot(figsize=(20,5))

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
ax = world.plot(figsize=(20,5))
ax.axis('off')

world.crs

ax = world.plot(figsize=(20,5), linewidth=0.25, edgecolor='white', color='lightgrey')
ax.axis('off')

ax = districts.to_crs(epsg=3857).plot(figsize=(20,5), linewidth=0.25, edgecolor='white', color='pink')
ax.axis('off')

earthquakes.plot(figsize=(20,7))

# Plot the earthquakes:
# BIGGER CIRCLES
# all one color
# make them a little transparent
# NO AXES OR BORDERS AROUND THE MAP!!!!!!!

ax = earthquakes.plot(figsize=(20,5), markersize=10, color='pink', alpha=0.5)
ax.axis('off')

# Plot the earthquakes:
# BIGGER CIRCLES
# all one color
# make them a little transparent
# NO AXES OR BORDERS AROUND THE MAP!!!!!!!

# Save the first layer as ax
ax = world.plot(color='lightgrey', linewidth=0.5, edgecolor='white', figsize=(15,5))
# Pass ax=ax to the second layer
earthquakes.plot(markersize=10, color='pink', alpha=0.5, ax=ax)
ax.axis('off')

earthquakes.head(2)

# Plot the earthquakes:
# BIGGER CIRCLES
# all one color
# make them a little transparent
# NO AXES OR BORDERS AROUND THE MAP!!!!!!!

# Save the first layer as ax
ax = world.plot(color='lightgrey', linewidth=0.5, edgecolor='white', figsize=(15,5))
# Pass ax=ax to the second layer
earthquakes.plot(markersize=10, alpha=0.5, ax=ax, column='mag')
ax.axis('off')

# Plot the earthquakes:
# BIGGER CIRCLES
# all one color
# make them a little transparent
# NO AXES OR BORDERS AROUND THE MAP!!!!!!!

# Save the first layer as ax
ax = world.plot(color='lightgrey', linewidth=0.5, edgecolor='white', figsize=(15,5))
# Pass ax=ax to the second layer
earthquakes.plot(markersize=10, alpha=0.5, ax=ax, column='mag', cmap='plasma')
ax.axis('off')

world.head(2)

world.plot(column='gdp_md_est', cmap='inferno')

world.plot(column='pop_est', cmap='inferno')

ax = world.plot()
ax.set_title("Default")

ax = world.to_crs({'proj': 'merc'}).plot()
ax.set_title("Default")

ax = world.to_crs({'proj': 'tmerc'}).plot()
ax.set_title("Default")

ax = world.to_crs({'proj': 'aea'}).plot()
ax.set_title("Default")

# Read in the shapefile from cb_2016_us_state_500k as "states"
states = gpd.read_file("cb_2016_us_state_500k/cb_2016_us_state_500k.shp")

# Get rid of Guam, Mariana Islands and Virgin Islands
states = states[states.STATEFP.astype(int) < 60]
# Get rid of Hawaii and Alaska
states = states[~states.NAME.isin(['Hawaii', 'Alaska'])]
states.tail(5)

# Read in the CSV
df = pd.read_csv("../wafflehouses.csv")

points = df.apply(lambda row: Point(row.long, row.lat), axis=1)
wafflehouses = gpd.GeoDataFrame(df, geometry=points)
wafflehouses.crs = {'init': 'epsg:4326'}

# If you want to know how this all works, look above
wafflehouses.head(2)

ax = states.plot()
wafflehouses.plot(column='score', ax=ax)

states.crs

wafflehouses.crs

# Convert CRS to match
converted_states = states.to_crs(wafflehouses.crs)
converted_states.head(2)

converted_states.crs

wafflehouses.crs

wafflehouses.head()

converted_states.head(2)

converted_states.plot()

wafflehouses.plot()

# http://geopandas.org/mergingdata.html
wafflehouses_with_state_data = gpd.sjoin(wafflehouses, converted_states, how='left', op='within')
wafflehouses_with_state_data.head()

wafflehouses_with_state_data['NAME'].value_counts()

# Need dropna because some wafflehouses are missing states
wafflehouses_with_state_data.dropna().plot(column='NAME', markersize=10, figsize=(20,5))

# http://geopandas.org/mergingdata.html
states_with_wafflehouse_data = gpd.sjoin(converted_states, wafflehouses, how='left', op='contains')
states_with_wafflehouse_data.head()

states_with_wafflehouse_data.shape

select_columns = states_with_wafflehouse_data[['NAME', 'geometry', 'score']]
select_columns.head()

states.head(2)

states.loc[0].geometry

# False = 0
# True = 1
# Count the number of states where the name contains "New"
states.NAME.str.contains("New").sum()

# Give me the first state!
state = states.loc[0]
# Look at the wafflehouses...
# are they inside of the state's geometry?
wafflehouses.within(state.geometry).sum()

# So for our first state, there were 147 inside of there

# You can use .contains
# counts the true ones
states['wafflehouse_count'] = states.apply(lambda state: wafflehouses.within(state.geometry).sum(), axis=1)
states.head()

states.plot(column='wafflehouse_count', scheme='quantiles')

states.plot(column='wafflehouse_count', scheme='equal_interval')

states.plot(column='wafflehouse_count', scheme='fisher_jenks', legend=True)

wafflehouses_with_state_data.to_file("wafflehouses.json", driver='GeoJSON')