<jupyter_start><jupyter_text>Geospatial Data in Python Vector DataWe are going to use [geopandas](https://geopandas.org/) to work with some vector data layers.  Geopandas is built on top of pandas, the main addition is its ability to handle **Geometries**.  For this, it relies on another package called [shapely](https://shapely.readthedocs.io/en/stable/manual.html)<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import geopandas as gpd


%matplotlib inline

Y = np.array([2,12,18,12,2])
X = np.array([12,12,15,18,18])

fig,ax=plt.subplots()
ax.scatter(X,Y)
ax.grid()
ax.set_title('Verticies')
from shapely.geometry import Polygon, LineString, Point<jupyter_output><empty_output><jupyter_text>Shapely Objects Point(s)<jupyter_code>from shapely.geometry import Point

point = Point([X[0],Y[0]])

point
from shapely.geometry import MultiPoint

mulit_point = MultiPoint([x for x in zip(X, Y)])

mulit_point<jupyter_output><empty_output><jupyter_text>Polygon(s)<jupyter_code>from shapely.geometry import Polygon

poly = Polygon([coord for coord in zip(X, Y)])

poly
from shapely.geometry import MultiPolygon

Mpol = []
for i in range(0,61,20):
    Mpol.append(Polygon([coord for coord in zip(X+i, Y+i)]))

print(Mpol)
Mpoly = MultiPolygon(Mpol)

Mpoly<jupyter_output>[<shapely.geometry.polygon.Polygon object at 0x0000020B4E2F5640>, <shapely.geometry.polygon.Polygon object at 0x0000020B4E31AB80>, <shapely.geometry.polygon.Polygon object at 0x0000020B4A8649D0>, <shapely.geometry.polygon.Polygon object at 0x0000020B4A864A60>]<jupyter_text>Line(s)<jupyter_code>from shapely.geometry import LineString

line = LineString([x for x in zip(X, Y)])

line
from shapely.geometry import MultiLineString

Mlin = []
for i in range(0,61,20):
    Mlin.append(LineString([coord for coord in zip(X+i, Y+i)]))

print(Mlin)
Mline = MultiLineString(Mlin)

Mline<jupyter_output>[<shapely.geometry.linestring.LineString object at 0x0000020B4A86E0A0>, <shapely.geometry.linestring.LineString object at 0x0000020B4F3CEBB0>, <shapely.geometry.linestring.LineString object at 0x0000020B4A7E8070>, <shapely.geometry.linestring.LineString object at 0x0000020B4F4AFE80>]<jupyter_text>Unzip Data- Below is a .zip file of census sub-divisions from Simply Analytics.- We need to unzip it and inspect the metadata<jupyter_code>import zipfile
Shape_file='SimplyAnalytics_Shapefiles_2021-11-18_02_24_45_8483b3a2bdad04d76a370657ccde2675'
with zipfile.ZipFile('data/'+Shape_file+'.zip', 'r') as zip_ref:
    zip_ref.extractall('data/Census/')
# print()
meta_data = pd.read_csv('data/Census/variable_names.txt',header=None,sep='#')
print(meta_data.values)<jupyter_output>[['VALUE0\t'
  ' Household Population by Aboriginal Identity | Household Population For Aboriginal Identity | Aboriginal Identity, 2021']
 ['VALUE1\t' ' Basics | Total Population, 2021']
 ['VALUE2\t'
  ' Household Population by Visible Minority Status | Household Population For Visible Minority | Visible Minority Total, 2021']]<jupyter_text>Reading the shapefileWe can rely on [geopandas](https://geopandas.org/en/stable/), a spatial extension for pandas.  We can use Geopandas to read, manipulate, and write geospatial data.<jupyter_code>import geopandas as gpd

# the .read_file() function reads shapefiles
BC_Sub_Div = gpd.read_file('data/Census/'+Shape_file+'.shp')

meta_data = pd.read_csv('data/Census/variable_names.txt',header=None,sep='#')
print('Vrriable Names: \n',meta_data.values)

BC_Sub_Div<jupyter_output>Vrriable Names: 
 [['VALUE0\t'
  ' Household Population by Aboriginal Identity | Household Population For Aboriginal Identity | Aboriginal Identity, 2021']
 ['VALUE1\t' ' Basics | Total Population, 2021']
 ['VALUE2\t'
  ' Household Population by Visible Minority Status | Household Population For Visible Minority | Visible Minority Total, 2021']]<jupyter_text>Editing Data<jupyter_code>## Note - This is the terminology used by the census
## It can be found in the variable_name file that comes with the download
BC_Sub_Div = BC_Sub_Div.rename(columns={
'VALUE0': 'Indigenous Identity',
'VALUE1': 'Population',
'VALUE2': 'Visible Minority'
})

BC_Sub_Div.head()<jupyter_output><empty_output><jupyter_text>Plotting the Data- Display the data and inspect the projection<jupyter_code>fig,ax=plt.subplots(figsize=(8,8))
BC_Sub_Div.plot(ax=ax)
ax.legend()

BC_Sub_Div.crs<jupyter_output>No handles with labels found to put in legend.<jupyter_text>Re-project the Data<jupyter_code>BC_Sub_Div_Albers = BC_Sub_Div.to_crs('EPSG:3005')

fig,ax=plt.subplots(figsize=(8,8))
BC_Sub_Div_Albers.plot(ax=ax)
ax.legend()

BC_Sub_Div_Albers.crs<jupyter_output>No handles with labels found to put in legend.<jupyter_text>Inspect the Geometry<jupyter_code>BC_Sub_Div_Albers.area<jupyter_output><empty_output><jupyter_text>Map by a Column<jupyter_code>fig,ax=plt.subplots(figsize=(8,8))
BC_Sub_Div_Albers.plot(column='Population',ax=ax,legend=True,scheme="quantiles",edgecolor='k')<jupyter_output><empty_output><jupyter_text>Clipping Data<jupyter_code># BC_Boundary.to_file('Data/Census/BC_Boudary_File.shp')
BC_Boundary = gpd.read_file('Data/Census/BC_Boudary_File.shp')

fig,ax = plt.subplots(1,2,figsize=(8,5))
BC_Sub_Div_Albers.plot(edgecolor='k',ax=ax[0])
BC_Boundary.plot(edgecolor='k',ax=ax[1])

ax[0].set_title('Input Layer')
ax[1].set_title('Clip Layer')


BC_Sub_Div_Clip = gpd.clip(BC_Sub_Div_Albers,BC_Boundary)

fig,ax = plt.subplots(figsize=(5,7))
BC_Sub_Div_Clip.plot(edgecolor='k',ax=ax)
ax.set_title('Final Result')<jupyter_output><empty_output><jupyter_text>Calculate Population Density<jupyter_code>BC_Sub_Div_Clip['Pop_Density'] = BC_Sub_Div_Clip['Population']/BC_Sub_Div_Clip.area*1e6
BC_Sub_Div_Clip['Pop_Density']=BC_Sub_Div_Clip['Pop_Density'].fillna(0)

fig,ax=plt.subplots(figsize=(8,8))
BC_Sub_Div_Clip.plot(column='Pop_Density',ax=ax,legend=True,scheme="User_Defined", 
         classification_kwds=dict(bins=[
             BC_Sub_Div_Clip['Pop_Density'].quantile(.25),
             BC_Sub_Div_Clip['Pop_Density'].quantile(.5),
             BC_Sub_Div_Clip['Pop_Density'].quantile(.75),
             BC_Sub_Div_Clip['Pop_Density'].max()]),
                    edgecolor='k')

ax.set_xlim(1e6,1.4e6)
ax.set_ylim(4e5,6e5)

# BC_Sub_Div_Clip.to_file('data//BC_Sub_Divisions.shp')<jupyter_output><empty_output><jupyter_text>Data Classification<jupyter_code>BC_Sub_Div_Clip.loc[((BC_Sub_Div_Clip['Pop_Density']<=400)|(BC_Sub_Div_Clip['Population']<=1000)),
                    'Community_Type']='Rural'

BC_Sub_Div_Clip.loc[((BC_Sub_Div_Clip['Community_Type']!='Rural')&(BC_Sub_Div_Clip['Population']<30000)),
                    'Community_Type']='Small Poulation Center'

BC_Sub_Div_Clip.loc[((BC_Sub_Div_Clip['Community_Type']!='Rural')&(BC_Sub_Div_Clip['Population']>=100000)),
                    'Community_Type']='Large Urban Center'

BC_Sub_Div_Clip.loc[((BC_Sub_Div_Clip['Community_Type']!='Rural')&(BC_Sub_Div_Clip['Population']>=30000)&
                     (BC_Sub_Div_Clip['Population']<100000)),
                    'Community_Type']='Medium Poulation Center'

print('Community Type Summary BC')
print(BC_Sub_Div_Clip.groupby('Community_Type').count()['name'].sort_values())


fig,ax=plt.subplots(figsize=(8,8))
BC_Sub_Div_Clip.plot(column='Community_Type',ax=ax,legend=True,cmap='Pastel2',edgecolor='k',linewidth=.25)

ax.set_xlim(1e6,1.4e6)
ax.set_ylim(3.5e5,6e5)

ax.set_title('Community Type SW BC')<jupyter_output>Community Type Summary BC
Community_Type
Medium Poulation Center      9
Large Urban Center          11
Small Poulation Center      44
Rural                      673
Name: name, dtype: int64<jupyter_text>Saving Data ShapefilesWe can save new layer as a shapefile.  Note the **"UserWarning"*** [Shapefiles](https://en.wikipedia.org/wiki/Shapefile) are only one type of vector file.  They have certain restrictions to minimize storage space (eg. limiting column names), they aren't human readable, and they split the data into multiple files.<jupyter_code>BC_Sub_Div_Clip.to_file('data/Outputs/BC_Sub_Divisions.shp')<jupyter_output>C:\Users\User\AppData\Local\Temp/ipykernel_10560/165744027.py:1: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  BC_Sub_Div_Clip.to_file('data/Outputs/BC_Sub_Divisions.shp')<jupyter_text>GeojsonWe can save the new layer as a [geojson](https://en.wikipedia.org/wiki/GeoJSON) file to get around the issue.* This is a simple file type often used in web-centered applications that stores data as a **human readable** dictionary.* The file takes up more space, but is also a bit more accessible/flexible.<jupyter_code>BC_Sub_Div_Clip.to_file("data/Outputs/BC_Sub_Div.json", driver = "GeoJSON")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12550097331959234, 'fraction_numerical': 0.035726554448643075, 'mean_word_length': 6.634615384615385, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 67, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '643636', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3452, 'n_tokens_neox': 3179, 'n_words': 592, 'path': '_site/Python_Notebooks/Working with Geopandas.ipynb', 'repo_name': 'ubc-library-rc/Geospatial-Analysis-Visualization-with-Python', 'size': '8733'}