<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')<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 0x000002332E7E36D0>, <shapely.geometry.polygon.Polygon object at 0x000002332E7E3370>, <shapely.geometry.polygon.Polygon object at 0x000002332E7E3700>, <shapely.geometry.polygon.Polygon object at 0x000002332E7E3730>]<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 0x000002332E7435E0>, <shapely.geometry.linestring.LineString object at 0x000002332E7E3FA0>, <shapely.geometry.linestring.LineString object at 0x000002332E7E3B20>, <shapely.geometry.linestring.LineString object at 0x000002332E01C5B0>]<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_04_29_59_93f600838bff00a6da2283b90dbf31c8'
with zipfile.ZipFile('data/'+Shape_file+'.zip', 'r') as zip_ref:
    zip_ref.extractall('data/Census/')<jupyter_output><empty_output><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.* We can open .txt files with pandas to view the metadata<jupyter_code>import geopandas as gpd

# the .read_file() function reads shapefiles
BC_FSA = 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_FSA<jupyter_output>Vrriable Names: 
 [['VALUE0\t' ' Basics | Total Population, 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_FSA = BC_FSA.rename(columns={
'VALUE0': 'Population',
})

BC_FSA.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_FSA.plot(ax=ax)
ax.legend()

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

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

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

fig,ax = plt.subplots(1,2,figsize=(8,5))
BC_FSA_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')<jupyter_output><empty_output><jupyter_text>ClipThis will take a little while to process<jupyter_code>BC_FSA_Clip = gpd.clip(BC_FSA_Albers,BC_Boundary)

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

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

ax.grid()

# ## Set Zoom & Turn Grid off
x = 1.225e6
y = 0.45e6
v = 8.5e4 
h = 8.5e4

ax.set_xlim(x-h,x+h)
ax.set_ylim(y-h,y+h)
ax.grid()

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

BC_FSA_Clip.loc[(BC_FSA_Clip['Community_Type']!='Rural'),
                    'Community_Type']='Urban'

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


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

ax.set_xlim(x-h,x+h)
ax.set_ylim(y-h,y+h)
# ax.grid()

ax.set_title('Community Type SW BC')<jupyter_output>Community Type Summary BC
Community_Type
Rural     86
Urban    103
Name: name, dtype: int64<jupyter_text>Saving Data ShapefilesWe can save new layer as a shapefile.  [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, are restricted to one geometry type (Polygon, line, or point), and they split the data into multiple files.Note the **"UserWarnings"**  * The first one isn't too serious, but it explains why the data came with a generic column header    *<jupyter_code>print(BC_FSA_Clip.groupby(BC_FSA_Clip['geometry'].type).count()['Population'])
BC_FSA_Clip.loc[BC_FSA_Clip['geometry'].type !='GeometryCollection'].to_file('data/Outputs/BC_FSA.shp')<jupyter_output>GeometryCollection      1
MultiPolygon           61
Polygon               127
Name: Population, dtype: int64<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>print(BC_FSA_Clip.groupby(BC_FSA_Clip['geometry'].type).count()['Population'])
BC_FSA_Clip.to_file("data/Outputs/BC_FSA.json", driver = "GeoJSON")<jupyter_output>GeometryCollection      1
MultiPolygon           61
Polygon               127
Name: Population, dtype: int64
{'fraction_non_alphanumeric': 0.1289117836445589, 'fraction_numerical': 0.034642627173213134, 'mean_word_length': 6.460134486071086, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 69, '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': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '43365', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3133, 'n_tokens_neox': 2909, 'n_words': 555, 'path': 'Python_Notebooks/Working with Geopandas.ipynb', 'repo_name': 'ubc-library-rc/Geospatial-Analysis-Visualization-with-Python', 'size': '7765'}