<jupyter_start><jupyter_text>Accessibility to schools - Sindh, Pakistan Inputs- Driving friction surface- Gridded population- School dataThis notebook cleans the school data, and re-uses inputs from the previous health script (urban population, and travel time) to calculate rural pop. within 30 min. drive of schoolsNearly identical to the KP notebook, but some adjustments had to be made to work with the school data for Sindh<jupyter_code>import sys, os
import rasterio

import pandas as pd
import geopandas as gpd
import numpy as np
import skimage.graph as graph

from rasterio.mask import mask
from rasterio import features
from shapely.geometry import box, Point, Polygon
from scipy.ndimage import generic_filter

# sys.path.append("../../../GOST_Urban")
# import src.UrbanRaster as urban

sys.path.append("/home/wb514197/Repos/INFRA_SAP")
sys.path.append('/home/wb514197/Repos/GOSTnets')
import infrasap.market_access as ma
import infrasap.rasterMisc as rMisc
from infrasap.misc import tPrint

import fiona
import matplotlib.pyplot as plt
%load_ext autoreload
%autoreload 2
%matplotlib inline
#User records
# destinations = "/home/public/Data/COUNTRY/PAK/HEALTH/pakistan_health_facilities/Pakistan_Health_Facilities.shp"
iso3 = "PAK"
data_dir = "/home/wb514197/data/PAK"
out_folder = os.path.join(data_dir, "output")

if not os.path.exists(out_folder):
    os.makedirs(out_folder)<jupyter_output><empty_output><jupyter_text>Process and clean school data<jupyter_code>data_dir
schools_dir = os.path.join(data_dir, "sindh")
fiona.drvsupport.supported_drivers['LIBKML'] = 'rw'
# fiona.drvsupport.supported_drivers['KML'] = 'rw'
os.chdir(os.path.join(data_dir, 'sindh_kml'))
for district in os.listdir(schools_dir):
    with open(os.path.join(schools_dir, district), 'rt', encoding="utf-8") as myfile:
        doc=myfile.read()
        link = doc[doc.find("CDATA[")+6:doc.find("]]></href>\n")]
        print(link)
schools_dir = os.path.join(data_dir, "sindh_kml")
kmz.columns
res_list[0].description.iloc[0]
res_list = []
for district in os.listdir(schools_dir):
    print(district)
    kmz = gpd.read_file(os.path.join(schools_dir, district))
#     for idx, row in kmz.iterrows():
#         table = pd.read_html(kmz.loc[idx].description)[1]
#         table = table.transpose()
#         table.columns = table.iloc[0]
#         table = table.drop(axis=0, index=0)
#         table.index = [idx]
#         if idx==0:
#             res = table.copy()
#         else:
#             res.loc[idx] = table.loc[idx]
    res_list.append(kmz)
res_df = pd.concat(res_list)
geoms = [Point(xy) for xy in zip(res_df.X.astype('float'), res_df.Y.astype('float'))]
crs = 'EPSG:4326'
schools = gpd.GeoDataFrame(res_df, crs=crs, geometry=geoms)
schools = res_df.copy()
schools.head()
schools.set_index("FID", inplace=True)
schools.to_file(os.path.join(data_dir, "schools_sindh.shp"), driver="ESRI Shapefile")
# DATA FROM PREVIOUS SCRIPT
admin1 = "/home/public/Data/GLOBAL/ADMIN/Admin1_Polys.shp"
admin2 = "/home/wb514197/data/PAK/pakistan_indicators.shp"

out_health = "/home/wb514197/data/PAK/HEALTH/"
out_pop_surface_std = os.path.join(out_health, "POP_2020_NEW_STD.tif")
out_travel_surface = os.path.join(out_health, "TRAVEL_SURFACE.tif")
out_pop_urban = os.path.join(out_health, "URBAN_POP.tif")
urban_raster = rasterio.open(out_pop_urban)
rMisc.standardizeInputRasters(urban_raster, travel_surf, os.path.join(out_health, "URBAN_POP_STD.tif"), data_type="C")
urban_raster_std = rasterio.open(os.path.join(out_health, "URBAN_POP_STD.tif"))
urban_data = urban_raster_std.read(1)
inG = rasterio.open(out_travel_surface)
pop_raster = rasterio.open(out_pop_surface_std)
inB2 = gpd.read_file(admin2)
inB2 = inB2.to_crs(inG.crs)
inB2 = inB2.loc[inB2.Year==2018].copy()
inB2.index.name = "FID"
inB2.reset_index(inplace=True)
inB2.Province.unique()
inB2_sel = inB2.loc[inB2.Province=="Sindh"].copy()<jupyter_output><empty_output><jupyter_text>Get rural pop<jupyter_code>urban_bool_ = rasterio.open(os.path.join(out_health, "URBAN.tif"))
urban_bool = urban_bool_.read(1)
urban_bool[urban_bool==1] = 2
urban_bool[urban_bool==0] = 1
urban_bool[urban_bool==2] = 0
rural_bool = urban_bool.copy()
pop_ = rasterio.open(os.path.join(out_health, "POP_2020_NEW.tif"))
pop = pop_.read(1)
pop[pop<0] = 0
pop_rural = pop*rural_bool
meta = pop_.meta.copy()
out_file = os.path.join(out_health, "RURAL_POP.tif")
with rasterio.open(out_file, 'w', **meta) as out:
    out.write_band(1, pop_rural)
rural_surf = rasterio.open(out_file)
rMisc.standardizeInputRasters(rural_surf, travel_surf, os.path.join(out_health, "RURAL_POP_STD.tif"), data_type="C")
rural_pop_ras = rasterio.open(os.path.join(out_health, "RURAL_POP_STD.tif"))
rural_pop_data = rural_pop_ras.read(1)
rural_pop_data.max()<jupyter_output><empty_output><jupyter_text>Filter schools<jupyter_code>schools.columns
schools.loc[schools.School_Level=="Higher Sec.", "School_Level"] = "Secondary"
schools.loc[schools.School_Level=="Elementary", "School_Level"] = "Primary"
schools.loc[schools.SchoolLeve=="Higher Secondary", "SchoolLeve"] = "Secondary"
schools.loc[schools.SchoolLeve=="High", "SchoolLeve"] = "Secondary"
schools.School_Gender.value_counts()
# primary = schools.loc[schools.SchoolLeve=="Primary"].copy()
# middle = schools.loc[schools.SchoolLeve=="Middle"].copy()
# secondary = schools.loc[(schools.SchoolLeve=="High") | (schools.SchoolLeve=="Higher Secondary")].copy()
thresh = 30
# create MCP object
inG_data = inG.read(1) * 1000
# Correct no data values
inG_data[inG_data < 0] = 99999999
mcp = graph.MCP_Geometric(inG_data)
# mcp = graph.MCP_Geometric(inG_data[0,:,:])
# facility_cells = ma.get_mcp_dests(inG, schools_filt)
# costs, traceback = mcp.find_costs(facility_cells) 
tt_raster = ma.calculate_travel_time(inG, mcp, schools_filt)[0]<jupyter_output><empty_output><jupyter_text>Combine population and travel time<jupyter_code>pop_data = pop_raster.read(1)
urban_data = urban_raster_std.read(1)
urban_data[urban_data==-0.0] = 0
pop_data[pop_data<0] = 0
cur_tt_d = (tt_raster < thresh).astype('int')
cur_pop = pop_data * cur_tt_d
from rasterstats import zonal_stats
schools_filt.columns
final = []
rural_pop_ras
rural_pop_data
schools = schools_original.copy()
schools.columns
for school_type in schools.School_Level.unique():
    gender= 'All'
    col = school_type+', '+gender
    print(f"\t {col}")
    schools_filt = schools.loc[schools.School_Level==school_type].copy()
    tt_raster = ma.calculate_travel_time(inG, mcp, schools_filt)[0]
    cur_tt_d = (tt_raster < thresh).astype('int')
    cur_pop = rural_pop_data * cur_tt_d
    res = pd.DataFrame(zonal_stats(inB2_sel, cur_pop, affine=rural_pop_ras.transform, stats='sum', nodata=pop_raster.nodata)).rename(columns = {'sum':col})
    final.append(res)
    for gender in schools_filt.School_Gender.unique():
        col = school_type+', '+gender
        schools_filt2 = schools_filt.loc[schools_filt.School_Gender==gender].copy()
        tt_raster = ma.calculate_travel_time(inG, mcp, schools_filt2)[0]
        cur_tt_d = (tt_raster < thresh).astype('int')
        cur_pop = rural_pop_data * cur_tt_d
        res = pd.DataFrame(zonal_stats(inB2_sel, cur_pop, affine=rural_pop_ras.transform, stats='sum', nodata=pop_raster.nodata)).rename(columns = {'sum':col})
        print(f"\t {col}")
        final.append(res)
final = pd.concat(final, axis=1)
zs_pop = pd.DataFrame(zonal_stats(inB2_sel, rural_pop_data, affine=pop_raster.transform, stats='sum', nodata=rural_pop_ras.nodata)).rename(columns = {'sum':'rural_pop'})
inB2_sel.columns
output = inB2_sel[['FID', 'Province', 'District']].copy()
output = output.reset_index(drop=True)
output = output.join(zs_pop).join(final)
out_folder
output.to_csv(os.path.join(out_folder, "Share of Rural Pop Within 30 Min Sindh.csv"), index=False)
os.getcwd()
os.listdir()<jupyter_output><empty_output><jupyter_text>of students<jupyter_code>gpd.sjoin
inB2_sel.loc[:,"District2"] = inB2_sel.loc[:,"District"]
schools_sj = gpd.sjoin(schools, inB2_sel[['Province', 'District2', 'geometry']], how='left')
schools_original = schools.copy()
schools = schools_sj.copy()
schools.loc[:,'Boys'] = schools.loc[:,'Boys'].astype('int')
schools.loc[:,'Girls'] = schools.loc[:,'Girls'].astype('int')
schools.loc[(schools.Boys>0) & (schools.Girls>0), ['SchoolGend', 'Boys', 'Girls']]
schools.loc[(schools.Boys>0) & (schools.Girls>0), 'SchoolGend'] = "Mixed"
schools[['School_Gender', 'School_Leve','Boys','Girls']]
students = schools.groupby(['District', 'SchoolLeve'])['Boys', 'Girls'].sum()
#District
students2 = schools.groupby(['District2', 'SchoolLeve', 'SchoolGend'])['Boys', 'Girls'].sum()
students2
students.to_csv(os.path.join(out_folder, "Students by Gender.csv"), index=True)
students2.to_csv(os.path.join(out_folder, "Students by Gender and Mixed2.csv"), index=True)<jupyter_output><empty_output><jupyter_text>of schools<jupyter_code>schools_n = schools.groupby(['District2', 'SchoolLeve', 'SchoolGend'])[['Scode']].count()
schools_n.to_csv(os.path.join(out_folder, "Number of Schools by Gender and Mixed2.csv"), index=True)<jupyter_output><empty_output><jupyter_text>Average TT<jupyter_code>inG_data.max()
inG.transform
inG_data.shape
zs_avg = pd.DataFrame(zonal_stats(inB2, inG_data, affine=inG.transform, stats='mean', nodata=inG_data.max())).rename(columns = {'mean':'average travel time'})
pop_data.min()
zs_sum_pop = pd.DataFrame(zonal_stats(inB2, pop_data, affine=pop_raster.transform, stats='sum', nodata=pop_raster.nodata)).rename(columns={'sum':'pop_sum'})
target_temp = inB2.join(zs_sum_pop)
type(target_temp)
from infrasap import aggregator
pop_raster
pop_sum_array = aggregator.rasterize_gdf(inD=target_temp, field='pop_sum', template=out_pop_surface_std)
pop_data = pop_raster.read(1)
pop_data[pop_data<0] = np.nan
plt.imshow(pop_data)
pop_data.min()
weights = pop_data/pop_sum_array
weights
pop_sum_array.shape
plt.imshow(weights)
zs_sum_weights = pd.DataFrame(zonal_stats(inB2, weights, affine=pop_raster.transform, stats='sum', nodata=np.nan))
inG_data = inG.read(1, masked=False) * 1000
# Correct no data values
# inG_data[inG_data < 0] = 99999999
inG_data.shape
inG_data[inG_data < 0] = np.nan
inG_data.max()
data_weighted = weights*inG_data
zs_weighted = pd.DataFrame(zonal_stats(inB2, data_weighted, affine=inG.transform, stats='sum', nodata=np.nan)).rename(columns={'sum':'pop weighted average travel time'})
zs_weighted.max()
tt_table = inB2[['Province', 'District']].join(zs_avg).join(zs_weighted)
tt_table.to_csv(os.path.join(out_folder, "Average Travel Time.csv"), index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14921416398409393, 'fraction_numerical': 0.015527362241999621, 'mean_word_length': 6.507462686567164, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '109891', 'license': 'MIT', 'n_tokens_mistral': 4322, 'n_tokens_neox': 4113, 'n_words': 653, 'path': 'Implementations/FY21/ACC_Pakistan - Accessibility Health and Schools/PAK_Accessibility_Schools_Sindh.ipynb', 'repo_name': 'worldbank/GOST_PublicGoods', 'size': '10562'}