+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ import pandas as pd +import geopandas as gpd +from sqlalchemy import create_engine +import numpy as np ++
+
+
+
+
+
+
+
+
+
+
+
+ Loading and exploring data¶
+
+
+
+
+
+
+
+
+
+
+
+
+ Loading data¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Pull in a layer from the OpenSGID +connection_url = "postgresql+psycopg2://agrc:agrc@opensgid.ugrc.utah.gov:5432/opensgid" +connection = create_engine(connection_url) +sql = "SELECT * FROM recreation.parks_local" +parks = gpd.read_postgis(sql, connection, geom_col="shape") ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ parks
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Dataframes and Series¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Get info about the dataframe +parks.info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Show a specific column +parks["name"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Get info about a specific column +parks["name"].info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Use .loc to get a specific row +parks.loc[0] ++
+
+
+
+
+
+
+
+
+
+
+
+ The Index¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ parks.index ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Create copy of dataframe indexed by park name +parks_by_name = parks.set_index("name") +parks_by_name ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Check if the index is unique +parks_by_name.index.is_unique ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: View duplicated index values +parks_by_name[parks_by_name.index.duplicated(keep=False)].sort_index() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can view the shape of the dataframe to see how many rows and columns it has +parks_by_name.shape ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Passing the dataframe to the len() function will return the number of rows +len(parks_by_name) ++
+
+
+
+
+
+
+
+
+
+
+
+ Selecting Data¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: .loc: selection by label +parks_by_name.loc["Flat Iron Mesa Park"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Get multiple dimensions with .loc +parks_by_name.loc["Flat Iron Mesa Park", "city"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: .iloc: selection by position +parks.iloc[0] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: same since we have a range index 0 to n +parks.loc[0] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: errors since now '0' is not a valid index label +parks_by_name.loc[0] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Multiple dimensions with slicing +parks.iloc[0, 3:] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: []: shorthand for selecting columns by label +parks["city"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: all rows and column "city" +parks.loc[:, "city"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Can use a list of column names to select multiple columns +parks[["name", "city"]] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Can also use a list of booleans to select rows +parks_head = parks.head() +parks_head ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Pass list of booleans representing the rows we want to keep +parks_head[[True, False, True, True, True]] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can use an equality operator to create a series of booleans whose index matches the dataframe +boolean_series = parks_head["city"] == "Logan" +boolean_series ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ parks_head[boolean_series] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Do it all in one step +parks_head[parks_head["city"] == "Logan"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ parks[parks["county"] == "KANE"] ++
+
+
+
+
+
+
+
+
+
+
+
+ Where clauses¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: You can do where-style queries with .query() as well +parks.query("county == 'KANE'") ++
+
+
+
+
+
+
+
+
+
+
+
+ Exploring the data¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: What are the unique cities represented? +parks["city"].unique() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: How many unique cities are there (including Nulls)? +parks["city"].nunique(dropna=False) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: How many parks are in each city (including Nulls)? +parks["city"].value_counts(dropna=False) ++
+
+
+
+
+
+
+
+
+
+
+
+ Cleaning, Analyzing, and Manipulating Data¶
+
+
+
+
+
+
+
+
+
+
+
+
+ Data Types and Missing Data¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Read in a new dataset from a zipped shapefile +opportunity_zones = gpd.read_file("data/ozs_shp_utah.zip") +opportunity_zones.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We've got some columns with numbers, let's check their data types ("dtype") +opportunity_zones.info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: let's try converting mfi_rat to a float +opportunity_zones["mfi_rat"] = opportunity_zones["mfi_rat"].astype(float) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Hmm, 'not available'. Let's see how many values are 'not available' +opportunity_zones["mfi_rat"].value_counts() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We'll replace 'not available' with NaN, which stands for "Not a Number" and is the standard way to represent missing values in pandas. +opportunity_zones["mfi_rat"] = opportunity_zones["mfi_rat"].replace("not available", np.nan) +opportunity_zones["mfi_rat"] = opportunity_zones["mfi_rat"].astype(float) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ opportunity_zones.info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: how many values got replaced with nan? +#: Because python treats True as 1 and False as 0, we can sum a boolean series to count the number of True values. +opportunity_zones["mfi_rat"].isna().sum() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: What about mfi? +opportunity_zones["mfi"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's try converting to float +opportunity_zones["mfi"].astype(float) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We'll use string ops to remove the dollar sign then convert to float +opportunity_zones["mfi"].str.replace("$", "").astype(float) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's replace "not available" across the whole dataframe with NaN +opportunity_zones.replace("not available", np.nan, inplace=True) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now try the conversion +opportunity_zones["mfi"].str.replace("$", "").astype(float) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: What happens if we convert to int instead of float? +opportunity_zones["mfi"].str.replace("$", "").astype(int) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: np.nan is a floating point value. Default int doesn't understand missing data. +#: We can either use floats or panda's nullable integer type +opportunity_zones["mfi"].str.replace("$", "").astype("Int64") #: Note the capital "I" in "Int64" ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: What's up with "1e+05"? Let's find it +opportunity_zones[opportunity_zones["mfi"].str.contains("1e", na=False)] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: The best we can do here is to replace "1e+05" with "100000" before converting to numeric +opportunity_zones["mfi"].str.replace("1e\\+05", "100000", regex=True).str.replace("$", "").astype("Int64") ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now we can do it all in one long chained operation +#: Chained operations work when each operation returns a new series (or dataframe) +opportunity_zones["mfi"] = opportunity_zones["mfi"].str.replace("1e\\+05", "100000", regex=True).str.replace("$", "").astype("Int64") ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's verify those columns look right +opportunity_zones.info() ++
+
+
+
+
+
+
+
+
+
+
+
+ Analyzing Rows: Vectorized Operations¶
+
+
+
+
+
+
+
+
+
+
+
+
+ The key to understanding pandas operations is to think in terms of operations that are applied to + every element in a row (or column), rather than extracting each element and passing it to the + operation one by one.
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Standard python: we're responsible for looping through a vector and calling a function on each element: +for element in [1, 2, 3]: + print(np.sqrt(element)) ++
+
+
+
+
+
+
+
+
+
+
+
+ 
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Vectorized operation: pass a whole vector (group of values) to a vectorized function: +np.sqrt([1, 2, 3]) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: pandas operation: a method on a series or dataframe +pd.Series([1, 2, 3]).transform(np.sqrt) ++
+
+
+
+
+
+
+
+
+
+
+
+ 
+
+
+
+
+
+
+
+
+
+
+
+ for loops in python are computationally expensive and require extra
+ resources to set up the iteration. In addition, the function has to be called multiple times, requiring
+ even more work behind the scenes for each call.
In contrast, vectorized operations are optimized to perform the same operation on multiple pieces
+ of data. In addition to avoiding the overhead from iteration and multiple function calls, the
+ processor has special logic and routines for parallelizing many operations. However, to use these it needs
+ to know the operation and the data type ahead of time, which it generally can't with python
+ for loops.
+
+
+
+
+
+
+
+
+
+
+
+ 
+
+
+
+
+
+
+
+
+
+
+
+ Think + about sending a set of data to an operation, not operating on data one piece at a time.¶ +
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's round the mfi column to the nearest thousand +opportunity_zones["mfi"].round(-3) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: As we saw when cleaning the data, string operations are vectorized as well. Let's capitalize the county names: +opportunity_zones["COUNTY_"].str.upper() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: The previous examples are transform operations that take in a series, operate on each element, and return a series of the same size. +#: If we want to use all the values in a series to compute a single value, we can use an aggregation operation instead. +opportunity_zones["mfi"].agg("mean") ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Many operations are available directly on the series +opportunity_zones["mfi"].mean() ++
+
+
+
+
+
+
+
+
+
+
+
+ .apply() and lambda functions: Applying any + function¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: .apply() is a generalized way of applying a function to a series or dataframe. +#: According to the docs, this is the preferred method (though some examples show looping to be faster) +def round_to_nearest_thousand(x): + return round(x, -3) +opportunity_zones["mfi"].apply(round_to_nearest_thousand) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: lambda functions are small single-line functions that can be defined inline without a name. They can be useful for simple operations that don't need to be reused. +#: The syntax is 'lambda <arguments>: <expression>' +#: When we use apply on a series, it passes each element of the series as the argument (ie, 'x') +opportunity_zones["mfi"].apply(lambda x: round(x, -3)) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: If we want to reference values from multiple columns, we can apply the function to the whole dataframe, specify axis=1 to apply the function to each row, and reference the columns by name within the function. +opportunity_zones.apply(lambda x: x["mfi"] * x["mfi_rat"], axis=1).value_counts() ++
+
+
+
+
+
+
+
+
+
+
+
+ Groupby: Performing the same analysis + on different subgroups¶
+Think of groupby as the analysis equivalent to symbolizing by unique values.
+.groupby splits a dataframe by the values of a column,
+ applies an operation on that each chunk's sub-frame, and then combines
+ the results into a data structure based on the type of operation performed.
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: What's the average mfi by county? +opportunity_zones.groupby("COUNTY_")["mfi"].mean() ++
+
+
+
+
+
+
+
+
+
+
+
+ 
+
+
+
+
+
+
+
+
+
+
+
+ Each groupby chunk is its own DataFrame, and any operation that can be done on a DataFrame can be done to
+ the chunk. .groupby() returns a groupby object that handles the iteration over the
+ DataFrames, and it also gives you access to the individual groups' DataFrames
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Get a groupby object and list the groups +#: .groups is a dict, so let's just view the keys to see the group names +grouped = opportunity_zones.groupby('COUNTY_') +grouped.groups.keys() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ type(grouped) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ dir(grouped) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Access an individual group's dataframe +grouped.get_group('Sevier County').head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: In addition to the built-in functions, we can use .apply() to pass arbitrary functions (like lambda functions) row-wise +opportunity_zones.groupby("COUNTY_").apply(lambda x: x["mfi"].mean(), include_groups=False) ++
+
+
+
+
+
+
+
+
+
+
+
+ Reshaping¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's look at the eligibility of tracts by MSA +by_msa = opportunity_zones.groupby("msa")["oz_elig"].value_counts() +by_msa ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Our result is a series with a multi-index +type(by_msa) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ by_msa.index ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's make it look prettier with unstack() +by_msa.unstack() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Resetting the index gives us a dataframe with columns that contain the groups +by_msa.reset_index() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can accomplish the same thing as unstack() by pivoting the reset-indexed dataframe +by_msa.reset_index().pivot(index="msa", columns="oz_elig", values="count") ++
+
+
+
+
+
+
+
+
+
+
+
+ Joining and Updating Data¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Read in our csv with additional data +additional_data = pd.read_csv("data/additional_data.csv") +additional_data.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: To join, we need a common column to be our join key. We'll use the tract GEOID +additional_data["GEOID_tract"].info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Are the tract values unique? We don't want a one-to-many join +additional_data["GEOID_tract"].is_unique ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's check the type and uniqueness of the tract GEOID column in the opportunity zones dataframe as well +opportunity_zones["GEOID"].info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ opportunity_zones["GEOID"].is_unique ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: it's unique, but we've got a type mismatch. Let's convert the additional data GEOID to a string so it matches the opportunity zones GEOID +additional_data["GEOID_tract"] = additional_data["GEOID_tract"].astype(str) +additional_data["GEOID_tract"].info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Since we know both sets of GEOIDs are unique, we can use sets to check if they match before we join +set(additional_data["GEOID_tract"]) == set(opportunity_zones["GEOID"]) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now let's perform the join. We'll only use a couple of columns from the additional data +merged = opportunity_zones.merge(additional_data[["GEOID_tract", "OZ Eligibility", "Median Age 2024"]], left_on="GEOID", right_on="GEOID_tract", how="left") +merged.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now let's update the original eligibility column with the new data and drop the GEOID_tract and OZ Eligibility columns. This has the effect of updating data and adding a new column. +merged["oz_elig"] = merged["OZ Eligibility"] +merged.drop(columns=["GEOID_tract", "OZ Eligibility"], inplace=True) +merged.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now check to see if the values are different +merged["oz_elig"].value_counts() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ opportunity_zones["oz_elig"].value_counts() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can also use .update() to update values in place without merging. This needs the indices to match, so we'll set the index to the GEOID and then update the eligibility column with the new data + +#: we'll make some copies to preserve the originals +oz_copy = opportunity_zones.copy() +additional_copy = additional_data.copy() + +#: set the indices +oz_copy.set_index("GEOID", inplace=True) +additional_copy.set_index("GEOID_tract", inplace=True) + +#: use .update() to update the eligibility column in the opportunity zones dataframe with the new data. +#: This will only update values where the indices match, which is why we set them earlier. +#: We also need to rename the column in the view of the additional data so that it matches. +oz_copy.update(additional_copy[["OZ Eligibility"]].rename(columns={"OZ Eligibility": "oz_elig"})) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ oz_copy["oz_elig"].value_counts() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ additional_copy.info() ++
+
+
+
+
+
+
+
+
+
+
+
+ Let's Get Spatial¶
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's go back to our parks from the beginning +parks.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: This is a geodataframe- a dataframe with a geometry column and special spatial operations. +#: Let's look at the series currently set as the active geometry +parks.geometry ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can also check the coordinate system +parks.crs ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We'll make a copy of the dataframe and reproject it to web mercator +parks_mercator = parks.to_crs(epsg=3857) +parks_mercator.geometry ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We're going to get a subset that we can plot easier +sandy_parks = parks[parks["city"] == "Sandy City"] ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ sandy_parks.plot() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can do all the basic geo operations like buffering, etc. Let's buffer the parks by 500 meters +buffered_parks = sandy_parks.copy() +buffered_parks["shape"] = buffered_parks.geometry.buffer(100) +buffered_parks.plot() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ # park_centroids = sandy_parks.geometry.centroid +# park_centroids.plot() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's answer a question: How many parks are within 1km of an aquifer recharge area? +#: First, get the aquifer recharge data from the SGID +sql = "SELECT * FROM geoscience.aquifer_recharge_discharge_areas" +aquifers = gpd.read_postgis(sql, connection, geom_col="shape") +aquifers.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: What are the values we need to select recharge areas? +aquifers["zone"].value_counts() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ recharge_areas = aquifers[aquifers["zone"].str.contains("recharge", case=False, na=False)] +recharge_areas.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ recharge_areas.plot() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: A geodataframe can have multiple geometry columns with only one active at a time. Let's create a new geometry column with the parks buffered by 1km and set it as the active geometry. +parks["shape_buffer"] = parks.geometry.buffer(1000) +parks.set_geometry("shape_buffer", inplace=True) +parks.plot() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now we can do a spatial join so that we know which parks are within the limit and can access the aquifer area data. +#: We'll do a left join with parks as left, so we have all the original park data, but the parks that don't intersect have null values for the aquifer data. +parks_with_aquifer_info = gpd.sjoin(parks, recharge_areas, how="left", predicate="intersects") +parks_with_aquifer_info.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ parks_with_aquifer_info["zone"].value_counts(dropna=False) #: dropna=False keeps the null values in the count so we can see how many parks don't intersect with recharge areas at all ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now we can subset down to just the parks that are within the area +parks_with_aquifer_info[parks_with_aquifer_info["zone"].notna()].plot() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: If we look at the list of columns above, we see we've got two geometry columns. Let's make a copy using our subset and investigate the current geometry. +parks_in_recharge_areas = parks_with_aquifer_info[parks_with_aquifer_info["zone"].notna()].copy() +parks_in_recharge_areas.geometry.name ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Change the geometry column and drop the other +parks_in_recharge_areas = parks_in_recharge_areas.set_geometry("shape").drop(columns=["shape_buffer", "index_right", "xid_right"]) +parks_in_recharge_areas.geometry.name ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Verify the crs +parks_in_recharge_areas.crs ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Thanks to the pyogrio module, we can write this to a file gdb. If the specified gdb doesn't exist, it will be created. If the layer already exists, it will be overwritten. +parks_in_recharge_areas.to_file("data/output.gdb", layer="parks_with_aquifer_info", driver="OpenFileGDB", engine="pyogrio") ++
+
+
+
+
+
+
+
+
+
+
+
+ Working with the arcgis Library¶
+We've been using geopandas, but Esri's arcgis library gives us the similar Spatially-Enabled DataFrames + along with lots of methods for interacting with AGOL/Portal.
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: First, we need the arcgis library +import arcgis ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Read in a layer from ArcGIS Online +item_id = "62a6b78d53a943db906de0abd1c9174f" #: Utah City and Town Locations ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: First we need a GIS object that represents an AGOL org. Since this is a public item, we can use an anonymous GIS object +gis = arcgis.gis.GIS() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We can access any items in our portal. +item = gis.content.get(item_id) +item ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Theres so much you can do with the ArcGIS API for Python, but we're gonna stick to working with Spatially-Enabled DataFrames. +item.layers ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: We could either access the layer from the list of layers in the item, or from the FeatureLayer.fromitem() method. +layer = arcgis.features.FeatureLayer.fromitem(item) +layer ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: and we can access the layer's data as an sedf via the as_df parameter of the .query() method. +city_points = layer.query(where="1=1", as_df=True).head() +city_points.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Just like geodataframes, we can check the geographic info of the sedf. +#: All the spatial operations of sedf's are under the .spatial accessor +city_points.spatial.geometry_type ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: check the spatial reference +city_points.spatial.sr ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: reproject to UTM 12N +city_points.spatial.project(26912) +city_points.spatial.sr ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: If you're in a python environment with arcpy, you can write to gdb. +#: Otherwise, you can write to shapefiles. +city_points.spatial.to_featureclass("data/output.gdb/city_points") ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: If we were logged into a Portal with write access, we could publish the SEDF as a feature layer +city_points.spatial.to_featurelayer("my awesome layer", gis, folder="best stuff") ++
+
+
+
+
+
+
+
+
+
+
+
+ Bring It All Together¶
+How many students are in public schools that are within opportunity zones vs schools not in OZs, split by + Elementary, Junior High/Middle School, and High School?
+
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Bring in the schools from OpenSGID +sql = "SELECT * FROM society.schools_prek_to_12" +schools = gpd.read_postgis(sql, connection, geom_col="shape") +schools.head() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ schools.info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Change from the strings "True" and "False" to boolean +schools["charterschool"] = schools["charterschool"].apply(lambda x: True if x == "True" else False) +schools["privateschool"] = schools["privateschool"].apply(lambda x: True if x == "True" else False) +schools["onlineschool"] = schools["onlineschool"].apply(lambda x: True if x == "True" else False) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's narrow down to public schools +public_schools = schools[~(schools["charterschool"]) & ~(schools["privateschool"]) & ~(schools["onlineschool"])].copy() +public_schools ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ oz_copy.info() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now we can intersect our schools with the opportunity zones to copy data from that table to the schools table +#: We'll use oz_copy since it's got the updated eligibility data +#: Let's make sure the coordinate systems match +public_schools.crs.name ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ oz_copy.crs.name ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: reproject oz to the schools crs +oz_copy.to_crs(public_schools.crs, inplace=True) ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Now we can do the spatial join +schools_and_zones = public_schools.sjoin(oz_copy[["geometry", "oz_elig"]], how="left") +schools_and_zones ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ totals = schools_and_zones.groupby("oz_elig").apply(lambda x: x.groupby("schoollevel")["totalk12"].sum(), include_groups=False) +totals ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Unstack our totals to make it purdy +totals.unstack() ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Let's flip it long for wide to see side-by-side comparisons +totals.unstack().T ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ #: Write it out to csv +totals.unstack().T.to_csv("data/students_in_zones.csv") ++
+
+
+
+
+
+
+
+ In [ ]:
+
+
+
+
+
+
+ ++