<jupyter_start><jupyter_text>Project 3 Codebook Topic: Spatial Analysis of Exposure to Respiratory Hazards and School Performance in Riverside, California.This codebook contains a step by step guide to download, retrieve, subset and visualize data pooled from different sources (EPA EJscreen, NCES, SEDA, ACS).Time Frame for Analysis1. EJScreen data - 2020 2. NCES/Districts - 2018/20193. ACS- 20184. SEDA- Standardized test scores administered in 3rd through 8 th grade in mathematics and Reading Language Arts (RLA) over the 2008-09 through 2017-18 school years.**Note:Some cells have  infront of code(remove hashtag to run the code in the cells)**  1: Download software<jupyter_code># Import software

import pandas as pd
import geopandas as gpd 
import matplotlib.pyplot as plt 
import contextily as ctx
import numpy as np
import quilt3
from geopandas_view import view 
import libpysal as lps
import seaborn as sns
import tobler as tob
plt.rcParams['figure.figsize'] = [20, 10]<jupyter_output><empty_output><jupyter_text>2: Retrieve and Adjust Data  2.1 Environmental Justice Screen Data From EPA<jupyter_code># Retrieve EPA EJ Screen data from UCR CGS Quilt Bucket

b = quilt3.Bucket("s3://spatial-ucr")
b.fetch("epa/ejscreen/ejscreen_2020.parquet", "./ejscreen_2020.parquet"), # Might be a good idea to get 2018
ej = pd.read_parquet('ejscreen_2020.parquet')
ej.head()
# Rename EJ column to state GEOID
ej.rename(columns = {'ID' : 'GEOID'}, inplace = True)
# ej.head()
ej.columns
# Download USA Census block groups from the 2018 ACS via Quilt
geoms = gpd.read_parquet('s3://spatial-ucr/census/acs/acs_2018_bg.parquet')
geoms.head()
# Merge EJ and ACS data
ej = geoms.merge(ej, on='GEOID')
# ej.head()
# Filter EJ Screen data so it only displays CA 
ca_ej = ej[ej.GEOID.str.startswith('06')]
ca_ej.head()
# Filter out EJ Index for Air toxics respiratory hazard index
D_RESP_2 = ca_ej.D_RESP_2 
ca_ej.D_RESP_2 = D_RESP_2.replace(to_replace= "None", value=np.nan).astype(float)
# Create variable for Riverside, Orange and LA County using the fips code

riv_ej = ca_ej[ca_ej.GEOID.str.startswith("06065")] # Riverside County EJSCREEN
# Checking Data for Orange County
#riv_ej.head()
riv_ej.shape
# New dataframe for respiratory hazards index, minority percent , low income percent 
new_riverside = riv_ej[['D_RESP_2','MINORPCT','LOWINCPCT', 'RESP','geometry']]
new_riverside.head()
#new_riverside.groupby(by ='MINORPCT').sum()
#Statistics
new_riverside.describe()<jupyter_output><empty_output><jupyter_text>3: Merge District Data w/ EJ Data<jupyter_code># Upload school districts data

districts = gpd.read_parquet('s3://spatial-ucr/nces/districts/school_districts_1819.parquet')
# Subset of the california schools from the US data
CA_dist = districts[districts.STATEFP == "06"]
CA_dist.head()
# Locate Riverside Unified School District (RUSD)- better to make sure its in California

rusd_CA = districts[districts['STATEFP'].str.lower().str.contains('06')]
rusd = rusd_CA[rusd_CA['NAME'].str.lower().str.contains('riverside unified')]
rusd.head()

# Now, let's overlay the EPA EJ data (that was combined with ACS data) on top of the RUSD shape

rusd_ej = gpd.overlay(riv_ej,rusd, how='intersection')<jupyter_output><empty_output><jupyter_text>4: Combine School Locations with School District/EJ Map<jupyter_code># Download NCES school location data 
schools = gpd.read_parquet("s3://spatial-ucr/nces/schools/schools_1819.parquet")
#schools.head()
# Download SEDA learning outcomes data
seda = pd.read_csv("https://stacks.stanford.edu/file/druid:db586ns4974/seda_school_pool_gcs_4.0.csv",converters={"sedasch":str})
seda.sedasch=seda.sedasch.str.rjust(12, "0")
#seda.head()
# Convert NCES data into a GeoDataFrame
school_geoms = schools[['NCESSCH','CNTY','NMCNTY', 'geometry']]
school_geoms.head()
# Merge SEDA and NCES data
seda_merge = seda.merge(school_geoms, left_on="sedasch", right_on= "NCESSCH")
#seda_merge.head()
# Convert merged NCES/SEDA data into a GeoDataFrame and plot it
seda_merge= gpd.GeoDataFrame(seda_merge)
#seda_merge.plot()
# Subset data to only locate schools in Riverside County 
riv_schools = seda_merge[seda_merge['CNTY']=='06065']
#riv_schools.plot()
# Subset school data to find schools in RUSD

rusd_schools = gpd.overlay(riv_schools, rusd, how='intersection')<jupyter_output><empty_output><jupyter_text>5: Voronoi Polygons<jupyter_code># Subset of rusd schools with only mean scores
rusd_pts= rusd_schools[['gcs_mn_avg_ol','geometry']]
rusd_pts.head()
rusd_pts.plot(column='gcs_mn_avg_ol', legend=True)
# Subset for EJ screen for RUSD only containing Respiratory index and geometry

EJ_RUSD= rusd_ej[['D_RESP_2','geometry']]
EJ_RUSD.plot(column='D_RESP_2', cmap='Greens', scheme='Quantiles', k=3,edgecolor='grey',
        legend=True)<jupyter_output><empty_output><jupyter_text>5.1 Overlays  Approach 1 Spatial Join<jupyter_code># 1 for 1 mapping
base  = EJ_RUSD.plot(column='D_RESP_2', legend=True)
rusd_pts.plot(color='red', ax=base)
rusd_sch= gpd.sjoin(rusd_pts, EJ_RUSD, how='left', op='within')
rusd_sch.head()
rusd_pts['D_RESP_2'] = rusd_sch.D_RESP_2
rusd_pts.head()
rusd_pts.crs<jupyter_output><empty_output><jupyter_text>Approach 2 Areal Interpolation<jupyter_code>x = rusd_pts.geometry.x
y = rusd_pts.geometry.y
cents = np.array([x,y]).T
cents
schools_vd, school_cents = lps.cg.voronoi_frames(cents)
base = schools_vd.plot()
rusd_pts.plot(ax=base, color='red')
EJ_RUSD.crs
base = EJ_RUSD.geometry.boundary.plot(edgecolor='green')
schools_vd.plot(ax=base)
rusd_pts.plot(ax=base, color='red')
schools_vd, school_cents = lps.cg.voronoi_frames(cents, clip = EJ_RUSD.unary_union)
base = EJ_RUSD.geometry.boundary.plot(edgecolor='green')
schools_vd.plot(ax=base)
rusd_pts.plot(ax=base, color='red')
base = EJ_RUSD.plot(column='D_RESP_2')
schools_vd.geometry.boundary.plot(ax=base, edgecolor='red')
rusd_pts.plot(ax=base, color='red')<jupyter_output><empty_output><jupyter_text>Estimate the escore for a school using areal interpolation<jupyter_code>RESP = tob.area_weighted.area_interpolate(source_df=EJ_RUSD,
                                           target_df=schools_vd,
                                           intensive_variables=['D_RESP_2'])
RESP 
base = RESP.plot(column='D_RESP_2')
rusd_pts.plot(ax=base, color='red')
x = EJ_RUSD.to_crs(epsg=4269)
gdf_proj = x.to_crs(EJ_RUSD.crs)
gdf_proj.head()
boundary_shape = cascaded_union(x.geometry)
coords = points_to_coords(gdf_proj.geometry)
# Calculate Voronoi Regions
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1175755734467568, 'fraction_numerical': 0.02065927388728543, 'mean_word_length': 6.211391018619934, 'pattern_counts': {'":': 1, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '583297', 'license': 'CC0-1.0', 'n_tokens_mistral': 2595, 'n_tokens_neox': 2431, 'n_words': 556, 'path': 'Voronoi_Test.ipynb', 'repo_name': 'WawNun/p280s21project3', 'size': '6583'}