```{python} import sys print(sys.executable) print(sys.version) ``` # Urban Greenspace and Asthma Prevalence ## Get Started with Big Data Pipelines Vegetation can provide many ecosystem services in urban areas, such as cleaner air and water and flood mitigation. However, the results are mixed on relationships between a simple measurement of vegetation cover (such as average NDVI, a measurement of vegetation health) and human health. We do, however, find relationships between landscape metrics that attempt to quantify the connectivity and structure of greenspace and human health. These types of metrics include mean patch size, edge density, and fragmentation. Read More: Study by [Tsai et al. 2019](https://doi.org/10.1007/s11252-018-0813-3) on the relationship between edge density and life expectancy in Baltimore, MD. The authors also discuss the influence of scale (e.g. the resolution of the imagery) on these relationships, which is important for this case study. In this notebook, you will write code to calculate patch, edge, and fragmentation statistics about urban greenspace in Chicago. These statistics should be reflective of the connectivity and spread of urban greenspace, which are important for ecosystem function and access. You will then use a linear model to identify statistically significant correlations between the distribution of greenspace and health data compiled by the US Center for Disease Control. ### Working with larger-than-memory (big) data For this project, we are going to split up the green space (NDVI) data by census tract, because this matches the human health data from the CDC. If we were interested in the average NDVI at this spatial scale, we could easily use a source of multispectral data like Landsat (30m) or even MODIS (250m) without a noticeable impact on our results. However, because we need to know more about the structure of green space within each tract, we need higher resolution data. For that, we will access the [National Agricultural Imagery Program (NAIP)](https://naip-usdaonline.hub.arcgis.com/) data, which is taken for the continental US every few years at 1m resolution. That’s enough to see individual trees and cars! The main purpose of the NAIP data is, as the name suggests, agriculture. However, it’s also a great source for urban areas where lots is happening on a very small scale. The NAIP data for the City of Chicago takes up about 20GB of space. This amount of data is likely to crash your kernel if you try to load it all in at once. It also would be inconvenient to store on your hard drive so that you can load it in a bit at a time for analysis. Even if your are using a computer that would be able to handle this amount of data, imagine if you were analysing the entire United States over multiple years! To help with this problem, you will use cloud-based tools to calculate your statistics instead of downloading rasters to your computer or container. You can crop the data entirely in the cloud, thereby saving on your memory and internet connection, using [rioxarray](https://corteva.github.io). ### Check your work with testing! This notebook does not have pre-built tests. You will need to write your own test code to make sure everything is working the way that you want. For many operations, this will be as simple as creating a plot to check that all your data lines up spatially the way you were expecting, or printing values as you go. However, if you don’t test as you go, you are likely to end up with intractable problems with the final code. ## STEP 1: Set up your analysis Import necessary packages - Create **reproducible file paths** for your project file structure. - To use cloud-optimized [GeoTiffs](https://www.earthdata.nasa.gov/about/esdis/esco/standards-practices/geotiff), we recommend some settings to make sure your code does not get stopped by a momentary connection lapse – see the code cell below. ```{python} %pip install -q -e .. ``` ```{python} from landmapyr.initial import robust_code, create_data_dir ``` ```{python} robust_code() data_dir = create_data_dir('chicago-greenspace') ``` ## STEP 2: Create a site map We use the [Center for Disease Control (CDC) Places](https://www.cdc.gov/places) dataset for human health data to compare with vegetation. CDC Places also provides some modified census tracts, clipped to the city boundary, to go along with the health data. We start by downloading the matching geographic data, and then select the City of Chicago. You can obtain urls for the U.S. Census Tract shapefiles from the [TIGER service](https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html). You’ll notice that these URLs use the [state FIPS](https://www.bls.gov/respondents/mwr/electronic-data-interchange/appendix-d-usps-state-abbreviations-and-fips-codes.htm), which you can get installing and using the [us](https://pypi.org/project/us/) package. - Download the Census tract Shapefile that goes along with CDC places - `tract_path = shp_tract_path(data_dir, place='chicago-tract')` - `chi_tract_gdf = download_census_tract(tract_path, placename='Chicago')` - Use a **conditional statement** to cache your download. - `if not os.path.exists(tract_path):` - Use a **row filter** to select only the census tracts in Chicago - `tract_gdf[tract_gdf.PlaceName == placename]` - Use a **spatial join** to select only the Census tracts that lie at least partially within the City of Chicago boundary. - **_Not sure what this means._** There is no need to cache the full dataset – stick with your pared down version containing only Chicago. ```{python} from landmapyr.cdcplaces import shp_tract_path, download_census_tract ``` ```{python} # Set up the census tract path tract_path = shp_tract_path(data_dir, 'chicago-tract') chi_tract_gdf = download_census_tract(tract_path, 'Chicago') ``` Code to save HV plot: ```{python} #| eval: false import hvplot.pandas from landmapyr.hv_plots import hvplot_tract_gdf chi_tract_hv = hvplot_tract_gdf(chi_tract_gdf) hvplot.save(chi_tract_hv, "chi_tract.html") ``` **NOW NEED TO incorporate this image** **Reflect and Respond** What do you notice about the City of Chicago from the coarse satellite image? Is green space evenly distributed? What can you learn about Chicago from websites, scientific papers, or other resources that might help explain what you see in the site map? **WRITE YOUR CITY OF CHICAGO DATA DESCRIPTION AND CITATION HERE** ### Download census tracts and select your urban area - Download the Census tract Shapefile for the state of Illinois (IL). - Use a **conditional statement** to cache the download ## STEP 3: Access Asthma and Urban Greenspaces Data ### Access human health data The U.S. Center for Disease Control (CDC) provides a number of health variables through their [Places Dataset](https://www.cdc.gov/places/) that might be correlated with urban greenspace. For this assignment, start with adult asthma. Try to limit the data as much as possible for download. Selecting the state and county is one way to do this. - Access Places data with an API, but as with many APIs it is easier to test out your search before building a URL. Navigate to the [Places Census Tract Data Portal](https://data.cdc.gov/500-Cities-Places/PLACES-Local-Data-for-Better-Health-Census-Tract-D) and search for the data you want. - The data portal will make an API call for you, but there is a simpler, easier to read and modify way to form an API call. Check out to the [socrata documentation](https://dev.socrata.com/docs/endpoints) to see how. You can also find some limited examples and a list of available parameters for this API on [CDC Places SODA Consumer API Documentation](https://dev.socrata.com/foundry/data.cdc.gov/cwsq-ngmh) - Once you have formed your query, you may notice that you have exactly 1000 rows. The Places SODA API limits you to 1000 records in a download. Either narrow your search or check out the `&$limit=` parameter to increase the number of rows downloaded. You can find more information on the [Paging page of the SODA API documentation](https://dev.socrata.com/docs/paging) - You should also clean up this data by renaming the `'data_value'` to something descriptive, and possibly selecting a subset of columns. ```{python} from landmapyr.cdcplaces import download_cdc_disease, join_tract_cdc from landmapyr.plots import plot_gdfs_map from landmapyr.naip import naip_path, download_naip_scenes, ndvi_naip_df ``` ```{python} # Preview asthma data cdc_df = download_cdc_disease(data_dir, 'asthma') cdc_df ``` ### Join health data with census tract boundaries - Join the census tract `GeoDataFrame` with the asthma prevalence `DataFrame` using the `.merge()` method. - You will need to change the data type of one of the merge columns to match, e.g. using `.astype('int64')` - There are a few census tracts in Chicago that do not have data. You should be able to confirm that they are not listed through the interactive Places Data Portal. However, if you have large chunks of the city missing, it may mean that you need to expand the record limit for your download. ```{python} #| label: fig-healthmap #| fig-cap: "Asthma Disease Severity" tract_cdc_gdf = join_tract_cdc(chi_tract_gdf, cdc_df) plot_gdfs_map(tract_cdc_gdf, column='asthma') ``` **Reflect and Respond** Write a description and citation for the asthma prevalence data. Do you notice anything about the spatial distribution of asthma in Chicago? From your research on the city, what might be some potential causes of any patterns you see? **ADD YOUR CDC PLACES DESCRIPTION AND CITATION HERE** ### Get NAIP Data URLs NAIP data are freely available through the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/) [SpatioTemporal Access Catalog (STAC)](https://stacspec.org/en). - Access planetary computer catalog via `pystac_client.Client.open()` - Loop across census tract to `.search()` for NAP data - `collections=["naip"]` - `intersects=shapely.to_geojson(tract_geometry)` - `datetime=f"{year}"` - Access URL using `search.assets['image'].href` - Accumulate the urls in a `pd.DataFrame` or `dict` for later - Occasionally you may find that the STAC service is momentarily unavailable. You should include code that will retry the request up to 5 times when you get the `pystac_client.exceptions.APIError` > **Warning** > > As always – DO NOT try to write this loop all at once! Stick with one > step at a time, making sure to test your work. You also probably want > to add a `break` into your loop to stop the loop after a single > iteration. This will help prevent long waits during debugging. Download NAIP scenes if not done already. Might want special case if some index values already downloaded. ```{python} naip_index_path = naip_path(data_dir, 'chicago') %store -r chi_scenes_df try: chi_scenes_df except NameError: chi_scenes_df = download_naip_scenes(naip_index_path, tract_cdc_gdf) %store chi_scenes_df ``` ### Compute NDVI Index Statistics Calculate some metrics to get at different aspects of the distribution of greenspace in each census tract. These **fragmentation statistics** can be implemented with the [scipy](https://scipy.org/) package. Some examples include: - Percentage vegetation: percent pixels that exceed a vegetation threshold (.12 is common with Landsat) - Patches: average size of **patches**, or contiguous area, exceeding the vegetation threshold. Patches can be identified with the `label` function from `scipy.ndimage` - Edges: proportion of edge pixels among vegetated pixels. Edges can be identified by **convolving** the image with a **kernel** designed to detect pixels that are different from their surroundings. > **What is convolution?** > > Referring to differential equations, > convolution is an approximation of the > [Laplace transform](https://mathworld.wolfram.com/LaplaceTransform.html). > For the purposes of calculating edge density, convolution means that > we are taking all the possible 3x3 chunks for our image, and > multiplying it by the kernel: > > $$ > \text{Kernel} = > \begin{bmatrix} > 1 & 1 & 1 \\ > 1 & -8 & 1 \\ > 1 & 1 & 1 > \end{bmatrix} > $$ > > The result is a matrix the same size as the original, minus the > outermost edge. > If the center pixel is the same as the surroundings, > its value in the final matrix will be > $-8 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 = 0$. > If it is higher than the surroundings, the result will be negative, > and if it is lower than the surroundings, the result will be positive. > As such, the edge pixels of our patches will be negative. - Select census row(s) - Select a single row from the census tract `GeoDataFrame` using e.g. `.loc[[0]]`. - Select all the rows from your URL `DataFrame` that match the census tract. - For each URL, crop, merge, clip, and compute NDVI for that census tract. - Set a threshold to get a binary mask of vegetation. - Use sample code to compute the fragmentation statistics. Feel free to add any other statistics you think are relevant, but make sure to include a fraction vegetation, mean patch size, and edge density. ### Repeat for all tracts - Using a loop, for each Census Tract: - Using a loop, for each data URL: - Use `rioxarray` to open up a connection to the STAC asset, just like you would a file on your computer - Crop and then clip your data to the census tract boundary > HINT: check out the `.clip_box` parameter `auto_expand` and the `clip` parameter `all_touched` to make sure you don’t end up with an empty array - Compute NDVI for the tract - Merge the NDVI rasters - Compute: - total number of pixels within the tract - fraction of pixels with an NDVI greater than .12 within the tract (and any other statistics you would like to look at) - Accumulate the statistics in a file for later - Using a conditional statement, ensure that you do not run this computation if you have already saved values. You do not want to run this step many times, or have to restart from scratch! There are many approaches to this, but we actually recommend implementing your caching in the previous cell when you generate your dataframe of URLs, since that step can take a few minutes as well. However, the important thing to cache is the computation. Basically, URLs are saved via StoreMagic in `chi_scenes_df`, while NDVI computations are stored in CSV file at `naip_index_path`; after computations, `ndvi_naip_df()` (**need better name**) returns the `DataFrame` `ndvi_index_df`. ```{python} ndvi_index_df = ndvi_naip_df(naip_index_path, tract_cdc_gdf, chi_scenes_df) ``` ## STEP 4: Explore your data with plots ### Chloropleth plots Before running any statistical models on your data, check that your download worked. You should see differences in both median income and mean NDVI across the city. Create a plot that contains: - 2 side-by-side Chloropleth plots - Asthma prevelence on one and mean NDVI on the other - Make sure to include a title and labeled color bars ```{python} from landmapyr.naip import merge_ndvi_cdc from landmapyr.explore import var_trans, train_test from landmapyr.plots import plot_gdfs_map, plot_matrix, plot_train_test ``` ```{python} #| label: fig-gdfsmap #| fig-cap: Asthma and Edge Density ndvi_cdc_gdf = merge_ndvi_cdc(tract_cdc_gdf, ndvi_index_df) plot_gdfs_map(ndvi_cdc_gdf) ``` Do you see any similarities in your plots? Do you think there is a relationship between adult asthma and any of your vegetation statistics in Chicago? Relate your visualization to the research you have done (the context of your analysis) if applicable. **ADD YOUR PLOT DESCRIPTION HERE** ## STEP 5: Explore a linear ordinary least-squares regression ### Model description One way to find if there is a statistically significant relationship between asthma prevalence and greenspace metrics is to run a linear ordinary least squares (OLS) regression and measure how well it is able to predict asthma given your chosen fragmentation statistics. Before fitting an OLS regression, you should check that your data are appropriate for the model. Write a model description for the linear ordinary least-squares regression that touches on: - What assumptions are made about the data? - What is the objective of this model? - What metrics could you use to evaluate the fit? - Advantages and potential problems with choosing this model. **ADD YOUR CDC PLACES DESCRIPTION AND CITATION HERE** ### Data preparation When fitting statistical models, you should make sure that your data meet the model assumptions through a process of selection and/or transformation. You can select data in various ways: - Eliminate observations (rows) or variables (columns) with missing data - Select model matching how variables are correlated - linear models are not good at modeling circles - Select variables that explain most of variability in response You can transform data: - transform while preserving order to better follow a bell shape (normal) or stabilize variance - `log` to manage right skew [sklearn.preprocessing.PowerTransformer](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PowerTransformer.html) - automated transformation choice schemes: [Box-Cox](https://www.geeksforgeeks.org/box-cox-transformation-using-python/) and [Bickel-Doksum](https://doi.org/10.2307/2287831) - offset and/or rescale to overcome negative numbers or effects caused by different ranges (e.g. `log(x - min(x) + 1)` or `(x - min(x))/(max(x) - min(x))`) - dimension reduction (e.g. principle component analysis (PCA)) to reduce multicollinearity among the predictors > **Tip** > > Keep in mind that data transforms like a log transform or a PCA must > be reversed (or rescaled) after modeling for the results to be meaningful. - Use `hvplot.scatter_matrix()` to create an exploratory plot of data. - Adjust data to address issues with linear model assumptions. - Check for `NaN` values, and drop rows and/or columns. Use the `.dropna()` method to drop rows with `NaN` values. - Explain any data transformations or selections you made and why ```{python} #| label: fig-plotmatrix #| fig-cap: "Model Matrix" logndvi_cdc_gdf = var_trans(ndvi_cdc_gdf) plot_matrix(logndvi_cdc_gdf) ``` **EXPLAIN YOUR SELECTION AND TRANSFORMATION PROCESS HERE** ### Fit and Predict The [scikitlearn](https://scikit-learn.org/) library has a slightly different approach than many software packages, emphasizing generic model performance measures like cross-validation and importance over coefficient p-values and correlation. The `scikitlearn` approach generalizes more smoothly to machine learning (ML) models where the statistical significance is harder to derive mathematically. - Use the `scikitlearn` documentation or ChatGPT as a starting point, split your data into training and testing datasets with [train_test_split()](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html). - Fit a linear regression to your training data. - Use your fitted model to predict the testing values. - Plot the predicted values against the measured values. You can use the following plotting code as a starting point. ```{python} #| label: fig-traintest #| fig-cap: "Observed vs Predicted Asthma" logndvi_cdc_test, reg, logndvi_cdc_gdf = train_test(logndvi_cdc_gdf) plot_train_test(logndvi_cdc_test) ``` ### Spatial bias We always need to think about bias, or systematic error, in model results. Every model is going to have some error, but we’d like to see that error evenly distributed. When the error is systematic, it can be an indication that we are missing something important in the model. In geographic data, we generally expect adjacent places to be more similar than distant places ([spatial autocorrelation](https://mgimond.github.io/Spatial/spatial-autocorrelation.html)). Linear models can address location through trends (as part of expected mean) or autocorrelation (as part of the variance-covariance structure). - Compute the model error (`predicted - measured`) for all census tracts - Plot the error as a chloropleth map with a diverging color scheme - Looking at both of your error plots, what do you notice? What are some possible explanations for any bias you see in your model? ```{python} #| label: fig-resid plot_gdfs_map(logndvi_cdc_gdf, column=['asthma','resid','edge_density'], color=['Blues','RdBu','Greens']) ``` GeoViews code not shown: ```{python} #| eval: false import holoviews as hv from landmapyr.gvplot import gvplot_ndvi_index, gvplot_resid model_fit = gvplot_ndvi_index(ndvi_cdc_gdf) resid = gvplot_resid(logndvi_cdc_gdf, reg, yvar='asthma') models_gv = (model_fit[0] + resid + model_fit[1]) hv.save(models_gv, 'bigdata_model.html') ``` **Reflect and Respond** What do you notice about your model from looking at the error plots? What additional data, transformations, or model type might help improve your results? ## References Tsai, Wei-Lun, Yu-Fai Leung, Melissa R. McHale, Myron F. Floyd, and Brian J. Reich. 2019. “Relationships Between Urban Green Land Cover and Human Health at Different Spatial Resolutions.” Urban Ecosystems 22 (2): 315–24. . Bickel PJ, Doksum K (1981) "An Analysis of Transformation Revisited" JASA 76. . **DESCRIBE AND INTERPRET YOUR IMAGE HERE**