Philadelphia - Comparison of the Redlining Map and Urban Green Space

Overview

In this case study of Philadelphia, I will download satellite-based multispectral data for the City of Philadelphia, and compare that to the redlining map and results from the U.S. Census American Community Survey. This is to help answer the question of how might have racism or other societal forces influenced the distribution of urban green space in Philadelphia.

Philadelphia was chosen as the site of choice because I went to Bryn Mawr College for undergrad, which is just west of Philadelphia proper on the 'Main Line'. I spent a lot of time in Philadelphia as well as studying the city from an architectural, historical, sociological, and anthropological perspectives for my courses. There are racial and economic disparities in (and surrounding) the city as well as other persisting inequalities that can be seen in the built environment as well as social factors.

Philadelphia, like many other cities, has uneven access to urban green space that is rooted in racial discrimination and other injustices prior to the City Beautiful Movement and New Deal era.

Data Description

There were two sources of data used for this project - Redlining and Harmonized Landsat Sentinel L30 data. Redlining data is available as georectified rasters and spatial data which can be downloaded as a geopackage here. The HLS data is satellite based multispectral data in the form of a raster image. Together using the two data sources, I can evaluate the connection between the two and see if there is a relationship between redlining and the NDVI (Normalized Difference Vegetation Index). Use of these Spectral Indices and other methods will be elaborated on further in the Methods Description next.

• Redlining - Georectified Rasters and Spatial Data

The scans of redlining maps and area descriptions are from the City Survey Files, 1935-1940 which are publicly accessible. Below is the scan of the Redlining Map of Philadelphia for reference but the digitally accessible version can be found here. It should be noted that current City of Philadelphia has areas not included on this map, and the redlining map also includes areas that are outside of the current City of Philadelphia. The georectified raters and spatial data for redlining is the spatialization of historic documents made by the Home Owners' Loan Corporation which was discussed in the Site Description above.

The spatial data used in this case study was compiled by a team of collaborators at Mapping Inequality: Redlining in New Deal America which includes professionals and experts in the field, professors, and students - for a full list of collaborators on this 3rd version of the spatial data used in this case study please see [here] (https://dsl.richmond.edu/panorama/redlining/about) and the data citations at the end of this section. The georectified rasters and spatial data are available under a CC-BY-NC license. Philadelphia data will be selected from this geopackage for redlining data, many other U.S. cities are also available, but the data is limited depending on availability like if a redlining map was ever created for that city or area, and if that redlining map is able to be converted spatially. For example, Cincinnati, Ohio had a redlining map created but the original formal map cannot be found, so Mapping Inequality does not have spatially available data for this city.

This data is essential for studying data justice and how the effects of redlining are still present.

Redlining map from Philadelphia, PA courtesy of Mapping Inequality ( Nelson and Winling (2023) )

• NASA Harmonized Landsat Sentinal - Multispectral Data

HLS L30 (Harmonized Landsat and Sentinel-2) data in short, is the combined measurement from the NASA/USGS Landsat satellites (8 and 9) and the European Sentinels (2A and 2B). The combination of the two is what makes them 'harmonized' and this also means there's more frequent observations (images) of the land being taken (every 2-3 days) rather than less frequent if used separately (Jeffrey G. Masek, Junchang Ju). The user guide can be found here. The L30 part is in reference to the images being at a special 30-meter spatial resolution, you can read more about it here. These satellites have multiple sensors, each sensor is for a different reflective spectral band in the electromagnetic wavelength. A band alone cannot convey much, but the relationship between two or more bands does - these are normalized spectral indices such as NDBI, NDVI, and others. Information on more in detail about remote sensing can be found here. Bands layered can also output a True Color Image that uses the visible bands (red, green, and blue) which is what we typically associate with a photo we take with a camera or our phone, and a CIR (Color Infrared) Image (a.k.a. False Color Image) which uses different bands (Near Infrared - 'NIR', red, and green). The use of NIR here enhances the visualization of healthy vegetation which the NIR strongly reflects, making healthy vegetation appear bright red - basically enhances what the human eye cannot see.

To access the data NASA Worldview was used to search for Philadelphia specifically looking at the HLS L30 dataset. A free account will need to be created to access this data and be able to use it in your codespace. The module earthaccess and function search_data using parameters are used to search for the HLS L30 data according to the parameters set. Because this data comes from 'professional' sources (NASA, USGS, and European agencies) and gathered through satellites, there should be a degree of trust in the data, but it depends how you are using it and what you are using it for. I am interested in the NDVI to tell a story of urban green space as it relates to the redlining map of Philadelphia.

The data can vary depending on the day or days used. For example, certain days have more or complete cloud coverage and it is essentially useless to for this case study to pick that single day. It is best to get a 'clear' image or images where there is little to no cloud coverage for the best results. I had to search for potential clear days that looked to be promising date to use as the 'temporal' parameter when trying to access this data (will be seen later in this portfolio post). There will also be differences in the image based on time of year depending on what you are looking at. If interested in vegetation for the U.S., it's important to know the climate and environment because you would want a time of year where the vegetation is healthiest and would have the greatest % reflection; however, it depends on what you are using this data to accomplish. In this case study a general vibrant time of year for vegetation like summer time in Philadelphia is ideal because the reflectance for vegetation at urban green spaces would be high at this time in theory. Therefore, it would be more clear where these green spaces are and help in the comparison to the redlining areas.

The data provided on NASA Worldview can be used for many things beyond what this case study aims to do, and it is worth exploring on your own.

Data Citations:

Redlining:

Redlining map of Philadelphia, PA courtesy of Mapping Inequality

• https://dsl.richmond.edu/panorama/redlining/about

Redlining - georectified rasters and spatial data:

• Nelson, Robert K, and LaDale Winling. 2023. “Mapping Inequality: Redlining in New Deal America.” In American Panorama: An Atlas of United States History, edited by Robert K Nelson and Edward L. Ayers. https://dsl.richmond.edu/panorama/redlining.

HLSL30 Multispectral Data:

NASA HLS L30 Product Description

• https://hls.gsfc.nasa.gov/products-description/l30/

USGS HLS L30 V002 Description

• Masek,Jeffrey G. and Ju, Junchang. "HLSL30 v002: HLS Operational Land Imager Surface Reflectance and TOA Brightness Daily Global 30m." https://lpdaac.usgs.gov/products/hlsl30v002/

NASA Worldview - Multispectral Data

• https://worldview.earthdata.nasa.gov

Methods Description

Going off of the 'Principle #2' of Data Feminism is to Challenge Power from authors Catherine D’Ignazio and Lauren Klein in Data Feminism, in order to analyze and expose oppression, certain methods can be used. One way to challenge power, rather than secure it, is to understand history, culture, and context (D’Ignazio and Klein, 2020). My goal is to "focus on data justice, rather than data ethics alone, (which) can help to ensure that past inequities are not distilled into black-boxed algorithms that, like the redlining maps of the twentieth century, determine the course of people’s lives in the twenty-first" (D’Ignazio and Klein, 2020). The 'Site Description' section, will be much longer than typically done for this reason. Putting Philadelphia in context of the larger (racist) utopian ideals that led to the decisions made during the the City Beautiful Movement and the New Deal Era is necessary to promote data justice.

Going into the coding side of things - after downloading needed packages, creating a data directory for the project, then downloading the redlining data, the multispectral data needs to be processed. This data needs to be processed because all the spectral bands and all the spatial areas for Philadelphia are needed to do further analysis.

The first step to processing the multispectral data is to download the data through earthaccess, search, and open all the rasters, in this case 2 cover Philadelphia. Then, a function is created to load these rasters, crop them (to reduce the amount of data needed dramatically), and scale them (redlining data is vector, and multispectral is raster so they need to be converted to the same CRS and encode the no data values correctly). Then, a second function is created to process the cloud mask, which is to exclude clouds, adjacent to clouds, and cloud shadow, as well as water from the image, then apply the cloud mask to the opened fmask layer of the raster data, the redlining gdf and the bits to mask using this function.

Next, data is prepped to be processed repeatedly using a regular expression, and then processed. Part of this processing includes creating functions with for loops to create 'DRY' code that repeats itself on each tile and on each band rather than creating repetitive code that would need to be done for each and every band and spatial area, which would not only be long, but also time consuming and unnecessary. The functions and use of for loops could also be used on if there are plural dates, however for the purpose of this case study, using one singular date will still show the effects of redlining on people's lives in the 21st century, but using plural dates should be thought about depending on what is being researched. If I were to further this case study, I would want to look at different time periods and plural dates as it relates to shifts in either the adding or taking away of urban green spaces in Philadelphia, if those dates or time frames are available as multispectral data (HLS L30).

After the for loop, the tiles need to be merged and this is to make the next images, this includes a 'true color' image and a CIR or 'false color' image. While the true color image looks like a photo with a camera, the false color image while may not look like a typical photo, it is easier to visually see the vegetation in this type of image because of the high contrast. For this case study the CIR image is very helpful to see the urban green spaces and add to the analysis.

To go beyond just an image and create more depth in the analysis, the multispectral data (HLS L30) will be used to look at a normalized spectral index (which is a relationship of two or more bands), specifically NDVI (Normalized Difference Vegetation Index, (nir -red) - (nir + red) ). There are many kinds of normalized spectral indices, but NDVI was chosen because that reflects where the urban green spaces are to be able to run further analysis on the relationship of urban green spaces and redlining. After converting both the NDVI data array and the redlining geodataframe into the same CRS (Mercator), the NDVI can be overlayed with the redlining outline. This will be in the form of a plot and an interactive plot, the interactive plot is particularly helpful in that hovering over it shows the coordinates, so, in analysis rough coordinates could be given for the viewer to better interpret and understand a location being discussed.

However, while the visual of the overlay is helpful, calculating zonal stats can provide even further strength in analysis by summarizing NDVI across same geographic areas as the redlining info. Zonal stats are calculated by using the zonal_stats function after converting vector data to a raster mask through regionmask package. Then, plotting the regional statistics in linked subplots will show the possible connection between redlining and urban green space.

Fitting a model and evaluating it, is another way to see if redlining is related to NDVI by trying to correctly predict the redlining grade from the mean NDVI value. Further explanation of this is in the Model Description and Evaluation near the end.

Together, the multiple images and plots provide multiple ways to see the relationship between redlining and urban green spaces, some with more depth than others, but they build upon one another and reinforce each other.

Methods Citations:

• D’Ignazio, Catherine, and Klein, Lauren. 2020. “2. Collect, Analyze, Imagine, Teach.” In Data Feminism.

Site Description

Setting the Stage - Architecture and Urban Planning History and Theory

History is a series of reactions, whether that be to a the previous artistic movement, rejection of some idea, or reaction to current events. Knowing the context around what lead to the Redlining Maps created during the New Deal Era, gives meaning and further understanding to these events as a reaction to theories and ideas laid in urban planning.

Lewis Mumford (1895-1990), a scholar in the history of cities and urban planning practice, describes a city as "a theatre of social action" to which the design of the city is secondary to cities' relationship to the natural environment and human culture (Mumford, 91-95). Mumford is one of many theorists and scholars who emphasized the natural environment and green space in relationship to 'the city' during the City Beautiful Movement. One of his many influences was Ebenezer Howard (1850-1928), an Englishman who grew up during the Urban Revolution and was influenced by utopian ideals to create the 'Garden City' (Howard, 328-335). The Garden City was a concentric ring diagram of a theorized city where there is a "central park containing important public buildings and surrounded by a 'Crystal Palace' ring of retail stores ... encircled by a permanent agricultural greenbelt of some 5,000 acres" (Howard, 328). This was a protest against urban overcrowding that Howard witnessed during the Urban Revolution. Another inspiration to both Mumford and Howard, was Fredrick Law Olmstead (1822-1903), a landscape architect who laid out Central Park in New York City, Bryn Mawr College's campus, and scores of many other parks and even suburbs (Olmstead, 321-327). He was a champion for funding public parks and can be considered a pioneer for a new way of looking at the urban built environment - the need to "address the central problem of the relationship between nature and the built environment" (Olmstead, 321).

There are other planners, architects, and philosophers who had a different reaction to the Garden City of Howard and Olmstead's ideas but still influenced the future 'green ideals' of cities. One of these is the 'Father of the Modernist Movement', Le Corbusier (1887-1965), whose reaction to the Garden City was creating 'The City of Tomorrow' (Le Corbusier, 336-344). Rather than spreading out, Le Corbusier idealized this City of Tomorrow as constructed vertically, therefore increasing the parks, gardens, an open spaces on the ground (Le Corbusier, 339-341). He saw green spaces as essential to balancing the man made urban development. Le Corbusier was not alone in this concept, others like Frank Lloyd Wright who theorized the 'Broadacre City', have also influenced architecture and urban planning to this day.

Olmstead, Howard, Le Corbusier, and Wright all influenced the City Beautiful Movement (1890's-1920's); however, Daniel Burnham is more frequently credited with the start of the City Beautiful Movement. At the 1893 Columbian Exposition with his and Olmstead's model titled a 'White City' featured "lagoons, fountains, promenades, statues, and neoclassical buildings (which) showcased how comprehensive planning might relieve urban problems" (Nepa, 2016). Burnham later developed a plan for Chicago in 1909 based on this (Kaiser and Godschalk, 403). The significance of the 1893 Columbian Exposition should not be understated; the "Fair’s Midway Plaisance contained “living zoos” of people from Dahomey, Java, Egypt, Turkey, and more" which further concretized the othering of foreigners, particularly people of color, to be less than human (Graff, 2022). This further contributed to structural racism which was then embodied in The City Beautiful Movement. While these neoclassical buildings and green spaces were emphasized to be 'for all' this excluded people of color who were viewed as less than human, and was to benefit white, wealthy people in America. Structural racism has deeper roots in colonialism and the City Beautiful Movement could be considered another way to continue control and exclude people of color.

Redlining in the New Deal Era

This sets the stage for the New Deal Era legislation which used the ideas and theories of the City Beautiful Movement. Mumford is one of the scholars who helped spread the ideas from this Movement to America during the New Deal Era (1932-1938). The New Deal was a series of programs and legislation by Franklin D. Roosevelt (FDR) the 32nd President of the U.S. whose goal was to address the effects of the Great Depression (1929-1941) which included a housing crisis where hundreds of thousands were foreclosed on or evicted due to joblessness from the economic crisis and inability to pay mortgages or rent. This resulted in mass homelessness and bank failures.

Part of these programs was the creation of a federal agency - Home Owners' Loan Corporation (HOLC) which was tasked with grading riskiness of American neighborhoods for banks and other lenders who made mortgages (About Redlining). For each of these neighborhoods a description of the residents was included, "if those residents were African Americans or, to a lesser extent, immigrants or Jews, HOLC deemed them a threat to the stability of home values and described their presence as an 'infiltration'." This practice was called Redlining - Redlining was legal and practiced for decades. It dramatically affected the relative wealth as well as the health of different racial groups in America (Mapping Inequality).

The practice of redlining in this movement and era used zoning ordinances and other policy (housing and insurance policy) to exclude low-income people who were predominately Black or of races and ethnicities other than White) and could only live in certain areas due to inability to afford, not being granted loans when or where a white man would be granted one, etc. This was deliberate action taken by numerous people including politicians, city planners, insurance people, people in real estate, etc. to purposefully disinvest in non-White neighborhoods and areas.

One of many effects of redlining was the disinvestment in urban centres, of which many had a predominant population of people of color and immigrants after the 'white flight' to the suburbs. As a reaction this degradation of urban centres, the 'New Urbanism' movement started in the 1980's and 1990's (Congress for the New Urbanism, 356-359). New Urbanism "views disinvestment in central cities, the spread of placeless sprawl, increasing separation by race and income, environmental deterioration, loss of agricultural lands and wilderness, and the erosion of society's built heritage as one interrelated community-building challenge" (Congress for the new Urbanism, 357). One of the principles they created was a range of parks distributed within neighborhoods and "conservation areas and open lands should be used to define and connect different neighborhoods and districts" (Congress for the New Urbanism, 358). This movement pulled ideas from the World Commission on Environment and Development (The Bruntland Commission), who made a call to action for sustainable development (The Bruntland Commission, 351-355).

Redlining is rooted in ideals and theories from urban planning history. The theories of the City Beautiful Movement and the legislation and programs from the New Deal Era were not revolutionary in being completely new ideas, they were borrowed from deeply entrenched structural racism carried out by urban planners, philosophers, and other historical figures (whom almost entirely are white men of wealth). Despite these more recent movements continued emphasis on the importance of green space in urban planning, the disproportionate effect of redlining on urban green spaces continues.

Physical Orientation of Philadelphia

The grid Philadelphia is based on has 2 axes, going (off axis) north to south is Broad St, and going (off axis) west to east is Market St. This center or zero axis is in what is known today as 'Center City' neighborhood which is excluded from the redlining map of Philadelphia, being labeled as 'central commercial' which is still true to this day (the Redlining map of Philadelphia is in the Data Description section). Center City is home to major tourist attractions like the Liberty Bell, City Hall, Philadelphia Art Museum, Ben Franklin Institute, and many other museums and civic buildings created during the City Beautiful Movement and the New Deal Era. To the northeast, east, and southeast, the Delaware River borders the city. To the west, the Schuylkill River cuts roughly north to south through the city. Many of the parks and urban green spaces surround or are near these and other waterfront areas.

Relevant Philadelphia History

Land Acknowledgement: First Peoples and Indigenous Peoples are the original stewards of the land taken before, during, and beyond colonialism in the Americas. Much of recorded history has prioritized Euro American perspectives and experiences that have misrepresented, excluded, and erased Black, Indigenous, and People of Color perspectives and experiences (History Colorado, Grounding Virtues). First People's and Indigenous People's stewardship and connection to this land must be honored and respected.

Please read the land acknowledgement above the map before continuing reading, it is imperative to naming and confronting racialized systems that dominate history and our current era.

The Commonwealth of Pennsylvania (colony) was founded by Englishman William Penn (1644-1718) in 1681 advocating for religious freedom as Penn was a Quaker (a religious sect), and Philadelphia founded in 1682 to be situated between the Delaware River and the Schuylkill River. Penn commissioned A Portraiture of the City of Philadelphia - 1863 by Thomas Holme (shown below) based on astronomical observations prior to extensive surveying and understanding of the physical geography which was particularly marshy along the Schuylkill (Brin, 2016). One of the most prominent features of this map as it relates to green space are the five public squares "Centre Square (later the site of City Hall), Northeast Square (today’s Franklin Square), Northwest Square (Logan Square), Southwest Square (Rittenhouse Square), and Southeast Square (Washington Square)" (Barry Lawrence Ruderman Map Collection, Stanford libraries).

Fast forward to the City Beautiful Movement era (1890's - 1920's), W.E.B. du Bois a preeminent intellectual of his generation (was awarded the first PhD degree from Harvard to an African American, founded the NAACP, and much, much more), was appointed to study "the Negro problem" of Philadelphia using a case study of the Seventh Ward, the city's black ghetto (du Bois, 110-116). At the time, "Philadelphia had the largest and oldest settlement of African Americans in the northern United States" with a black population that had nearly tripled between 1860-1890 (du Bois, 111). Throughout his book The Philadelphia Negro, he continually mentions color prejudice which keeps black people held back from being true citizens or even being considered human (du Bois, 113). All of what may be thought of as unrelated acts of discrimination individually such as: not being treated equally at work, not being able to get higher positions even if higher education attained, not allowed to enter new lines of work, having to pay more for housing that is worse than what white people have, being insulted or discriminated against in public, etc., are collective discrimination and prejudice towards one group to diminish them and make them less than human (du Bois, 113-115).

These discriminatory acts were further exacerbated by the housing crisis during the Great Depression, where the scarcity of housing caused Black enclaves to become further overcrowded, over-priced, and dilapidated (Taketomo, Mapping Inequality). This led to disastorous consequences - "December 1936, two rowhouses at 517-519 S 15th St in the 7th Ward collapsed, killing seven and injuring many more" (Taketomo, Mapping Inequality). To read the full description on what happened please visit Contexts: Black Philly After The Philadelphia Negro. The collapse of these rowhouses were catalysts for the City of Philadelphia, to raze (to completely destroy) whole blocks of the Seventh Ward using a $10 million dollar grant given from letter writing campaign by Black Philadelphians to the Federal Public Works Administration to get new housing (which didn't happen because the money, was used to destory housing) (Hunter, 2014). In 1937 the Philadelphia Housing Authority was created and received a $20 million dollar grant to create new housing, targeting two sites in North Philadelphia to be Black housing projects using the "racial composition rule" - this required "that housing projects be built in already heavily black areas, and that their residents should reflect the existing racial composition of the surrounding neighborhood" (Hunter, 2014). This displaced many Black Philadelphians to the north in the years to come, and today, what was then the Seventh Ward, a historically black neighborhood, almost nothing remains (Hunter, 2014). While there are historical markers (signs) in this area today, this was a deliberate use of housing policy to displace black residents and erasure of their roots and both physcical and social history in this neighborhood.

At the same time as the housing crisis (during the Great Depression) and the deliberate displacement of Black residents, the New Deal Era meant new parks and green spaces and civic buildings following the City Beautiful Movement. One of the major projects was Fairmount Parkway, later named Benjamin Franklin Parkway (1937), was made to create a ceremonial boulevard to Fairmount Park lined with civic institutions (Lynn). While some plans and discussions of it go back to the 1870's, there were many phases of plans, but it wasn't until 1907 when Paul Cret during the City Beautiful Movement, was influenced by Ecole des Beaux Arts Paris and the building of the "Avenue des Champs-Elysées and other great Parisian boulevards built during the reign of Napoleon III" to design Fairmount Parkway (Lynn). The parkway officially opened in 1918 and then between 1927 and 1941 many civic and arts buildings (Free Library, Rodin Museum, Philadelphia Art Museum, Municipal Court Building) (Lynn) were built along the parkway to impose grandiosity on the city and centralize these institutions into one area (Jacobs, 3-25).

The ideals of Olmstead, Howard, Mumford, Burnham, Le Corbusier, Wright and others that influenced the City Beautiful Movement and New Deal era, were taken as truth and not questioned. This was at the expense of people and communities, particularly Black communities and other communities of color or 'other-ness' (Jews, Immigrants, etc.), that saw slum clearance, mass displacement, and loss of physical cultural history.

Site Description Citations (in order by first mention):

• Mumford, Lewis, 1937. "What is a City?". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic,(pages 91-95).

• Howard, Ebenezer, 1898/1902. "'Author's Introduction' and 'The Town-Country Magnet' from 'Garden Cities of To-morrow'". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic,(pages 328-335).

• Olmstead, Fredrick Law, 1870. "'Public Parks and the Enlargement of Towns' from 'American Social Science Association'". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic, (pages 321-327).

• Le Corbusier, 1929. "'A Contemporary City' from 'The City of Tomorrow'". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic, (pages 336-344).

• Kaiser, Edward J. and Godschalk, David R., 1995. "'Twentieth-Century Land Use Planning and its Planning:A Stalwart Family Tree' from 'Journal of the American Planning Association'". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic, (pages 399-422).

• Nepa, Stephen, 2016 Rutgers University. "City Beautiful Movement". The Encyclopedia of Greater Philadelphia. https://philadelphiaencyclopedia.org/essays/city-beautiful-movement/

• Graff, Rebecca S., 2022. "The City Beyond the White City: The World's Columbian Exposition". Society of Architectural Historians & Charnley-Persky House Museum Foundation. (https://www.beyondthewhitecity.org/worlds-columbian-exposition#:~:text=The%20Fair's%20Midway%20Plaisance%20contained,and%20the%20Captive%20Balloon%20ride.)

• Congress for the New Urbanism, 1993. "Charter of the new Urbanism". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic, (pages 336-344).

• World Commission on Environment and Development (The Bruntland Commission), 1987. "'Towards Sustainable Development' from 'Our Common Future'". The City Reader, 5th Edition, edited by LeGates, T. Richard and Stout, Frederic, (pages 336-344).

• About Redlining - courtesy of Mapping Inequality (https://dsl.richmond.edu/panorama/redlining/about)

• History Colorado. "Grounding Virtues" from History Colorado Anti-Racism Work (https://www.historycolorado.org/sites/default/files/media/document/2020/Anti-Racism_Grounding_Virtues.pdf).

• Brin, Joseph G., 2015. "Courtyard Compass Reveals A City Off Axis". Hidden City: Exploring Philadelphia's Urban Landscape. (https://hiddencityphila.org/2015/08/courtyard-compass-reveals-a-city-off-axis/)

• Barry Lawrence Ruderman Map Collection, Stanford Libraries. Notes from "A Portraiture of the City of Philadelphia in the Province of Pennsylvania in America" with "A Letter from William Penn Proprietary and Governour of Pennsylvania in America, To The Committee of the Free Society of Traders" 1683. (https://exhibits.stanford.edu/ruderman/catalog/jj125yt6481).

• Hoffman, J.S., 2023. "Hotter, Wetter, Sneezier, & Wheezier: Present-day Environmental Disparity Among HOLC Neighborhoods." In Mapping Inequality: Redlining in New Deal America. (https://dsl.richmond.edu/panorama/redlining/environment)

• Taketomo, Kristian Mapping Inequality - Redlining in new Deal America - Philadelphia Context

• Hunter, Marcus Anthony, 2014. Contexts: "Black Philly After The Philadelphia Negro" from Contexts Winter 2014 Issue. (https://contexts.org/articles/black-philly-after-the-philadelphia-negro/).

• Miller, Lynn. "Benjamin Franklin Parkway" from the Encylopedia of Greater Philadelphia (https://philadelphiaencyclopedia.org/essays/benjamin-franklin-parkway/).

• Jacobs, Jane, 1961, Vintage Books Edition - Random House Inc. New York 1992. "The Death and Life of Great American Cities". (Pages 3-25).

Outline of project

A) Set up analysis

1. Import all packages needed
2. Define and create project data directory

B) Download redlining data

1. Define info for redlining download
2. Only download once
3. Load from file
4. Create interactive site map using hvplot

C) Download raster data through earthaccess, login, search, and load

D) Process raster image

1. Create function with description to process raster images
    a. Connect to the raster image
    b. Get the study bounds
    c. Crop
    d. Test the function

2. Create function with description to process the cloud mask
    a. Open fmask file
    b. Unpack the cloud mask bits
    c. Create list for bits to mask
    d. Define cloud_mask using np.sum
    e. Define philadelphia_cloud_mask using the function from 1 (above)
    f. Test the function

E) Prep for processing rater data repeatedly - Get metadata with regular expressions

1. Compile a regular expression to search for metadata
2. Find all the metadata in the file name
3. Create a DataFrame with the metadata
4. Add the File-like URI to the DataFrame

F) Process data using a for loop

1. Create dictionary with labels, for each band to process
2. Initialize structure for saving images - for loop
    a. Load cloud mask
    b. create another for loop, for the band_id, in the tile_df
        i. Create 'if' statement if band_id in band dictionary from (1) above
        to process and mask the band and then save in a data array for later
3. Merge all tiles

G) Plot true color and false color images

1. Plot a true color image
2. Plot a CIR image 

E) Plot raster (NDVI) with vector (redlining) overlay

1. Prepare the plot by reprojecting both NDVI and redlining gdf to 
the same CRS
2. Plot NDVI raster with redlining overlay using matplotlib
3. Plot NDVI raster with redlining overlay using hvplot

F) Summarize NDVI across same geographic areas the redlining info by caluculating zonal stats, plot to see possible connection between vegetation health and redlining

1. Convert vector data to a raster mask through regionmask package
2. Calculate zonal stats through the zonal_stats function
3. Plot the regional statistics in linked subplots. Remember to 
    a. Merge the NDVI stats with redlining geometry into one gdf
    b. Change grade to ordered Categorical for plotting
    c. Drop NA grades - like commercial areas that aren't given
    a redlining grade or category

G) Fit a Model

1. Fit a model - convert categories to numbers first
2. Plot model results
3. Evaluate the model using cross validation

# Set Up Analysis Part 1 of 2
## Import packages that will help with...

# Reproducible file paths
import os # Reproducible file paths
import pathlib # Find the home folder
import re # Extract metadata from file names
import zipfile # Work with zip files
from io import BytesIO # Stream binary (zip) files
from glob import glob # Find files by pattern

# Find files by pattern
import numpy as np # Unpack bit-wise Fmask - Process bit-wise cloud mask- adjust images 
import matplotlib.pyplot as plt # Overlay pandas and xarry plots,Overlay raster and vector data
import requests # Request data over HTTP
import rioxarray as rxr # Work with geospatial raster data

# Work with tabular, vector, and raster data
import cartopy.crs as ccrs # CRSs (Coordinate Reference Systems)
import geopandas as gpd # work with vector data
import hvplot.pandas # Interactive tabular and vector data
import hvplot.xarray # Interactive raster
import pandas as pd # Group and aggregate
import xarray as xr # Adjust images
from rioxarray.merge import merge_arrays # Merge rasters

# Access NASA data
import earthaccess # Access NASA data from the cloud

#  Work with zonal stats, decision tree, and evaluator for model
import regionmask # Convert shapefile to mask
from xrspatial import zonal_stats # Calculate zonal statistics

# Model to set threshholds in order to classify
from sklearn.tree import DecisionTreeClassifier, plot_tree # see what the threholds are
from sklearn.model_selection import train_test_split, cross_val_score # evaluate model

/opt/miniconda3/envs/earth-analytics-python/lib/python3.11/site-packages/dask/dataframe/__init__.py:42: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

# Set Up Analysis Part 2 of 2

# Define and create the project data directory
redlining_data_dir = os.path.join(
    pathlib.Path.home(),
    'earth-analytics',
    'data',
    'redlining'
)
os.makedirs(redlining_data_dir, exist_ok=True)

# Download Redlining Data Part 1 of 1

# Define info for redlining download
redlining_url = (
    "https://dsl.richmond.edu/panorama/redlining/static"
    "/mappinginequality.gpkg"
)
redlining_dir = os.path.join(redlining_data_dir, 'redlining')
os.makedirs(redlining_dir, exist_ok=True)
redlining_path = os.path.join(redlining_dir, 'redlining.shp')

# Only download once
if not os.path.exists(redlining_path):
    redlining_gdf = gpd.read_file(redlining_url)
    redlining_gdf.to_file(redlining_path)

# Load from file
redlining_gdf = gpd.read_file(redlining_path)

/opt/miniconda3/envs/earth-analytics-python/lib/python3.11/site-packages/pyogrio/raw.py:198: RuntimeWarning: /Users/briannagleason/earth-analytics/data/redlining/redlining/redlining.shp contains polygon(s) with rings with invalid winding order. Autocorrecting them, but that shapefile should be corrected using ogr2ogr for example.
  return ogr_read(

# Create an interactive site map, select data from Philadelphia
philadelphia_redlining_gdf = redlining_gdf[redlining_gdf.city=='Philadelphia']
philadelphia_redlining_gdf.dissolve().hvplot(
    geo=True, tiles='EsriImagery',
    title='City of Philadelphia - Redlining Outline',
    fill_color=None, line_color='darkorange', line_width=3,
    frame_width=600
)

Rough outlines of some urban green spaces can be seen outside of the central part of the city

The boxes or odd shapes within the boundary are urban green spaces or parks. These outlines show some of the major or flagship parks like Fairmount Park being the largest (follows the outline of the Schuylkill River, roughly -75.2 Longitude, 40 Latitude) and Penny Pack Park which follows the Penny Pack Creek ( roughly -75.1 Longitude, 40.025 Latitude). There are multiple 'squares' or smaller parks throughout the city such as: Rittenhouse Square, Logan Square, Franklin Square, and Washington Square, which are all in Center City but seem not to be visually captured in the outline. Some major ones are also not outlined because they fall outside of what areas were redlined for Philadelphia, but are part of the City of Philadelphia today (FDR Park which is a visibly green area just above the Delaware river but slightly northwest of the Schuylkill River, roughly -75.18 Longitude, 39.9 Latitude). This plot also conveys a similar non- contiguous areas as the redlining map does - not all areas are connected. Many of Philadelphia urban green spaces are centered around physical natural waterways like rivers, creeks, lakes, and ponds.

Some of these 'parks' are parks in that they are public spaces, but have no or little vegetation (sometimes as planters versus planted vegetation) and rather have water fountains or other water features. Some of these are 'Love Park' which is the iconic 'LOVE' sculpture wasn't installed until 1976; and Rail Park which was created this century takes use of old rail infrastructure to create public places. These 'parks' are relatively new and created post redlining and the City Beautiful and New Deal Eras.

This plot will help situate or orient the reader for the further analyses.

Citation: Visit Philly- Best Parks In Philadelphia

# Download Raster data through earthaccess Part 1 of 1

# login to earthaccess
earthaccess.login(strategy="interactive", persist=True)

# Search earthaccess
philadelphia_results = earthaccess.search_data(
    short_name="HLSL30",
    bounding_box=tuple(philadelphia_redlining_gdf.total_bounds),
    temporal=("2024-07-17", "2024-07-18")
)
philadelphia_results

# Open earthaccess results
philadelphia_files = earthaccess.open(philadelphia_results)

QUEUEING TASKS | :   0%|          | 0/30 [00:00<?, ?it/s]

PROCESSING TASKS | :   0%|          | 0/30 [00:00<?, ?it/s]

COLLECTING RESULTS | :   0%|          | 0/30 [00:00<?, ?it/s]

# Check philadelphia_files by calling it
## Will need to test a certain band in the next cell
## by using count of location in the list
philadelphia_files

[<File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.Fmask.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.VAA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.SZA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B10.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B01.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B04.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.VZA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B03.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.SAA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B06.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B07.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B09.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B02.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B11.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TVK.2024199T153924.v2.0/HLS.L30.T18TVK.2024199T153924.v2.0.B05.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.VZA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.SAA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B05.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B06.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.SZA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B10.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B07.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B04.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B01.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B11.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.Fmask.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B03.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B09.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.VAA.tif>,
 <File-like object HTTPFileSystem, https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T18TWK.2024199T153924.v2.0/HLS.L30.T18TWK.2024199T153924.v2.0.B02.tif>]

# Process Raster Image Part 1 of 2

# Create function with description to process raster images
def process_image(uri, bounds_gdf):
    """
    Load, crop, and scale a raster image from earthaccess

    Parameters
    ----------
    uri: file-like or path-like
      File accessor downloaded or obtained from earthaccess
    bounds_gdf: gpd.GeoDataFrame
      Area of interest to crop to

    Returns
    -------
    cropped_da: rxr.DataArray
      Processed raster
    """

    # Connect to the raster image
    da = rxr.open_rasterio(uri, mask_and_scale=True).squeeze()

    # Get the study bounds
    bounds = (
    bounds_gdf
    .to_crs(da.rio.crs)
    .total_bounds 
    )

    # Crop 
    cropped_da = (
    da.rio.clip_box(*bounds)
    )
    
    return cropped_da


# Test the function
process_image(philadelphia_files[12],philadelphia_redlining_gdf).plot(
    cbar_kwargs={"label": "Reflectance"},
    robust=True,
)
plt.gca().set(
    title='Philadelphia Processed Blue Band',
    xlabel='', 
    ylabel='', 
    xticks=[], 
    yticks=[]
)
plt.show()

# Process Raster Image Part 2 of 2

# Create function with description to process the cloud mask
def process_cloud_mask(cloud_uri, bounds_gdf, bits_to_mask):
    """
    Load an 8-bit Fmask file and process to a boolean mask

    Parameters
    ----------
    uri: file-like or path-like
      Fmask file accessor downloaded or obtained from earthaccess
    bounds_gdf: gpd.GeoDataFrame
      Area of interest to crop to
    bits_to_mask: list of int
      The indices of the bits to mask if set

    Returns
    -------
    cloud_mask: np.array
      Cloud mask
    """

    # Open fmask file 
    fmask_da = process_image(cloud_uri, bounds_gdf)

    # Unpack the cloud mask bits
    cloud_bits = (
    np.unpackbits(
        (
            # Get the cloud mask as an array...
            fmask_da.values
            # ... of 8-bit integers
            .astype('uint8')
            # With an extra axis to unpack the bits into
            [:, :, np.newaxis]
        ), 
        # List the least significant bit first to match the user guide
        bitorder='little',
        # Expand the array in a new dimension
        axis=-1)
    )
    
    # Create list for bits to mask
    bits_to_mask = [
    1, # Cloud
    2, # Adjacent to cloud
    3, # Cloud shadow
    5] # Water

    cloud_mask = np.sum(
    # Select bits 
    cloud_bits[:,:,bits_to_mask], 
    # Sum along the bit axis
    axis=-1
    # Check if any of bits 1, 2, or 3 are true
    ) == 0

    return cloud_mask


#  Define Philadelphia cloud mask
philadelphia_cloud_mask = process_cloud_mask(
    philadelphia_files[0], 
    philadelphia_redlining_gdf, 
    bits_to_mask= [1, 2, 3, 5] 
  )

# Test the function - chose a band and define that data array and plot
blue_da = process_image(philadelphia_files[12], philadelphia_redlining_gdf)
blue_da.where(philadelphia_cloud_mask).plot.imshow(
    cbar_kwargs={"label": "Reflectance"},
    robust=True,
)
plt.gca().set(
    title='Philadelphia Processed Blue Band with Cloud Mask',
    xlabel='', 
    ylabel='', 
    xticks=[], 
    yticks=[]
)
plt.show()

/var/folders/sy/dpb41wnn0w76w40gwfy01bv80000gn/T/ipykernel_92708/1155284852.py:33: RuntimeWarning: invalid value encountered in cast
  .astype('uint8')

# Prep for processing raster data repeatedly

# Compile a regular expression to search for metadata
metadata_rows = re.compile(
 r"HLS\.L30\.(?P<tile_id>T[0-9A-Z]+)\.(?P<date>\d+)T\d+\.v2\.0\."
 r"(?P<band_id>.+)\.tif"   
)

# Find all the metadata in the file name
philadelphia_metadata_rows = [
    metadata_rows.search(philadelphia_file.full_name).groupdict()
    for philadelphia_file in philadelphia_files]

# Create a DataFrame with the metadata
philadelphia_raster_df = pd.DataFrame(philadelphia_metadata_rows)
# Add the File-like URI to the DataFrame
philadelphia_raster_df['file'] = philadelphia_files
# Check the results
philadelphia_raster_df

	tile_id	date	band_id	file
0	T18TVK	2024199	Fmask	<File-like object HTTPFileSystem, https://data...
1	T18TVK	2024199	VAA	<File-like object HTTPFileSystem, https://data...
2	T18TVK	2024199	SZA	<File-like object HTTPFileSystem, https://data...
3	T18TVK	2024199	B10	<File-like object HTTPFileSystem, https://data...
4	T18TVK	2024199	B01	<File-like object HTTPFileSystem, https://data...
5	T18TVK	2024199	B04	<File-like object HTTPFileSystem, https://data...
6	T18TVK	2024199	VZA	<File-like object HTTPFileSystem, https://data...
7	T18TVK	2024199	B03	<File-like object HTTPFileSystem, https://data...
8	T18TVK	2024199	SAA	<File-like object HTTPFileSystem, https://data...
9	T18TVK	2024199	B06	<File-like object HTTPFileSystem, https://data...
10	T18TVK	2024199	B07	<File-like object HTTPFileSystem, https://data...
11	T18TVK	2024199	B09	<File-like object HTTPFileSystem, https://data...
12	T18TVK	2024199	B02	<File-like object HTTPFileSystem, https://data...
13	T18TVK	2024199	B11	<File-like object HTTPFileSystem, https://data...
14	T18TVK	2024199	B05	<File-like object HTTPFileSystem, https://data...
15	T18TWK	2024199	VZA	<File-like object HTTPFileSystem, https://data...
16	T18TWK	2024199	SAA	<File-like object HTTPFileSystem, https://data...
17	T18TWK	2024199	B05	<File-like object HTTPFileSystem, https://data...
18	T18TWK	2024199	B06	<File-like object HTTPFileSystem, https://data...
19	T18TWK	2024199	SZA	<File-like object HTTPFileSystem, https://data...
20	T18TWK	2024199	B10	<File-like object HTTPFileSystem, https://data...
21	T18TWK	2024199	B07	<File-like object HTTPFileSystem, https://data...
22	T18TWK	2024199	B04	<File-like object HTTPFileSystem, https://data...
23	T18TWK	2024199	B01	<File-like object HTTPFileSystem, https://data...
24	T18TWK	2024199	B11	<File-like object HTTPFileSystem, https://data...
25	T18TWK	2024199	Fmask	<File-like object HTTPFileSystem, https://data...
26	T18TWK	2024199	B03	<File-like object HTTPFileSystem, https://data...
27	T18TWK	2024199	B09	<File-like object HTTPFileSystem, https://data...
28	T18TWK	2024199	VAA	<File-like object HTTPFileSystem, https://data...
29	T18TWK	2024199	B02	<File-like object HTTPFileSystem, https://data...

# Process data using a for loop 

# Labels for each band to process - dictionary
philadelphia_bands = {
    'B02': 'blue',
    'B03': 'green',
    'B04': 'red',
    'B05': 'nir',
    'B06': 'swir1',
    'B07': 'swir2',
}
# Initialize structure for saving images
philadelphia_das = {band_name: [] for band_name in philadelphia_bands.values()}
for tile_id, tile_df in philadelphia_raster_df.groupby('tile_id'):

   # Load the cloud mask
   philadelphia_fmask_file = tile_df[tile_df.band_id=='Fmask'].file.values[0]
   philadelphia_cloud_mask = process_cloud_mask(
      philadelphia_fmask_file, 
      philadelphia_redlining_gdf, 
      [1, 2, 3, 5])

   for band_id, row in tile_df.groupby('band_id'):
      
        if band_id in philadelphia_bands:
            band_name = philadelphia_bands[band_id]
            # Process band
            philadelphia_band_da = process_image(
                row.file.values[0], 
                philadelphia_redlining_gdf)
           
            # Mask band
            philadelphia_band_masked_da = philadelphia_band_da.where(philadelphia_cloud_mask)

            # Store the resulting DataArray for later
            philadelphia_das[band_name].append(philadelphia_band_masked_da)  

/var/folders/sy/dpb41wnn0w76w40gwfy01bv80000gn/T/ipykernel_92708/1155284852.py:33: RuntimeWarning: invalid value encountered in cast
  .astype('uint8')

# Merge all tiles
philadelphia_merged_das = {
   band_name: merge_arrays(das) 
    for band_name, das 
    in philadelphia_das.items()}

# Plotting a True Color Image Part 1 of 1

# Define the rgb data array
philadelphia_rgb_da = (
    xr.concat(
        [
            philadelphia_merged_das['red'],
            philadelphia_merged_das['green'],
            philadelphia_merged_das['blue']
        ],
        dim='rgb')
)

# Use the stretch_rgb function with parameters and make descriptions
def stretch_rgb(rgb_da, low, high, brighten):
    """
    Short description: 
    Default images are too light/dark, stretch_rgb() can increase the visual contrast and brighten the image.

    Long description:
    Applies a linear stretch based on the low and high percentiles, which are used as the endpoints.

    Parameters
    ----------
    rgb_da: array-like
      Single da with red, green, and blue band with 'rgb' dimension defined prior to using function
    low: float or int
      Sets the low percentile endpoint for the stretch. The lower the percentile, the lighter the image.
    high: float or int
      Sets the high percentile endpoint for the stretch. The higher the percentile, the darker the image.
    brighten: float 
      Brightens the image, range is from 0-1 with 0 being no added brightness and 1 being full brightness.
  
    Returns
    -------
    rgb_da: array-like
      Outliers in brightness taken out of da -> 'streched' image
    """
    # Adjust low, high and brighten numbers as needed
    p_low, p_high = np.nanpercentile(rgb_da, (low, high))
    rgb_da = (rgb_da - p_low)  / (p_high - p_low) + brighten
    rgb_da = rgb_da.clip(0, 1)
    return rgb_da

philadelphia_rgb_da = stretch_rgb(philadelphia_rgb_da, 3, 96, .1)

philadelphia_rgb_da.hvplot.rgb(
    geo=True, 
    x='x', 
    y='y', 
    bands='rgb',
    xaxis=None, 
    yaxis=None,
    title='True Color (3 band) Philadelphia Image',
    width=600
 )

WARNING:param.GeoRGBPlot01801: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

True Color Image shows key waterways and larger green spaces in Philadelphia

This image, which has the cloud mask applied (black boxes in the lower left corner), shows the two key waterways that outline or create the City of Philadelphia's current border. The Delaware River which is the largest which cuts from southwest to northeast, and the Schuylkill River which is a feeder river to the Delaware cuts from the northwest to directly south; the current City of Philadelphia extends slightly southwest of this river similar to the Redlining Map. There are many other minor waterways that are not visible here. Some of the major green spaces like Fairmount Park (on both sides of the Schuylkill) and Penny Pack Park in the northeast. While this image is also more clear than the ESRI imagry tile from the outline of 'City of Philadelphia - Redlining Oultine' earlier in this post, however it is difficult to discern some of the smaller green spaces that are in the city.

# Plotting a CIR (False Color) Image Part 1 of 1

# Define the rgb data array
philadelphia_cir_da = (
    xr.concat(
        [
            philadelphia_merged_das['nir'] * .75,
            philadelphia_merged_das['red'],
            philadelphia_merged_das['green']
        ],
        dim='rgb')
)

philadelphia_cir_da = stretch_rgb(philadelphia_cir_da, 3, 96, .1)

#cir_da.plot.imshow(rgb='rgb')
philadelphia_cir_da.hvplot.rgb(
    geo=True, 
    x='x', 
    y='y', 
    bands='rgb',
    xaxis=None, 
    yaxis=None,
    title='CIR - False Color (3 band) Philadelphia Image',
    width=600
    )

WARNING:param.GeoRGBPlot01880: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

CIR (False Color) Image - smaller waterways and green spaces more visible than in True Color Image

Philadelphia while it is a very built up/dense city, there are many smaller green spaces that are easily seen in the CIR Image that are not seen or not as easily seen in the True Color Image. This is due to the enhanced visualization of healthy vegetation appearing red (use of the 'NIR' band). All of the major waterways and green spaces in this image are still present, but some of the smaller green spaces like the 'squares' in central city can be seen as well as very small green spaces throughout the city. This adds to the analysis that the larger green spaces are towards the larger waterways, which are at or near the edges of the city, and what green spaces are in the more central and dense parts of the city, are much smaller and dispersed.

Next, Redlining areas will be overlayed on NDVI for further analysis.

# Plot NDVI only - Part 1 of 1

# Create new da for NDVI then test plot

# Create new da using the Philadelphia_merged_das
philadelphia_ndvi_da = (
    (philadelphia_merged_das['nir'] - philadelphia_merged_das['red'])
    /
    (philadelphia_merged_das['nir'] + philadelphia_merged_das['red'])
)

# Plot NDVI to test
philadelphia_ndvi_da.plot(
    robust=True,
    cbar_kwargs={"label": " NDVI"}
    )
plt.gca().set(
    title='Philadelphia NDVI',
    xlabel='', 
    ylabel='', 
    xticks=[], 
    yticks=[]
)
plt.show()

# Plot raster (NDVI) with vector (redlining) overlay Part 1 of 3

# Prepare the Plot 
# Reproject Denver to Mercator
philadelphia_redlining_plot_gdf = philadelphia_redlining_gdf.to_crs(ccrs.Mercator())

# Reproject NDVI data and band raster data to Mercator
philadelphia_ndvi_plot_da = philadelphia_ndvi_da.rio.reproject(ccrs.Mercator())

# Plot raster (NDVI) with vector (redlining) overlay Part 2 of 3

# Plot raster with overlay xarray
plt.figure(figsize=(12, 8))

philadelphia_ndvi_plot_da.plot(
    vmin=0,
    robust=True,
    cmap='inferno',
    cbar_kwargs={"label": "NDVI"}
    )

philadelphia_redlining_plot_gdf.plot(
    ax=plt.gca(),
    color='none',
    edgecolor='white',
    linewidth=1.75
    ) 
                  
plt.gca().set(
    title='Philadelphia NDVI with Redlining Overlay', 
    xlabel='',
    ylabel='',
    xticks=[], 
    yticks=[],
    )
    
plt.show()

NDVI and Redlining Overlay (matplotlib) - NDVI highest near green spaces around waterways

The highest NDVI is a yellow or light yellow color on the above plot. Based on the site description and the previous images shown, Fairmount Park (center left) on both sides of the Schuylkill and Penny Pack Park (center, slightly northeast) are known green spaces in Philadelphia and also near or around waterways. Both of these green spaces have a very drastic difference NDVI compared to the areas around it that have much lower NDVI (light to dark purple on the plot). It is difficult in this plot to see smaller and fewer green spaces that are in the central part of the city, and what can be seen shows few orange to red areas, which are mid NDVI on the scale. Comparing to the Redlining map it appears these larger green spaces are physcial boundaries between better and worse graded Redlining areas. For example, to the north and west of Fairmount Park are the 'better' graded Redlining areas (grade/category A and B), versus the 'worse' grade Redlining areas (grade/ category C and D) to the south and east that have lower NDVI.

# Plot raster (NDVI) with vector (redlining) overlay Part 3 of 3

# Create interactive plot of raster with redlining overlay
(
    philadelphia_ndvi_plot_da.hvplot(
        robust=True,
        geo=True,
        xaxis=None, 
        yaxis=None,
        title='Interactive Philadelphia NDVI with Redlining Overlay',
        cmap='oranges',
        clabel='NDVI'
    )
    * philadelphia_redlining_plot_gdf.hvplot(
        geo=True,
        crs=ccrs.Mercator(),
        fill_color=None,
        line_width=2.5,
        width=500
        )
)

WARNING:param.GeoOverlayPlot02049: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

NDVI and Redlining Overlay (hvplot)- similar to matplotlib conclusion, NDVI highest near green spaces around waterways

Both matplotlib and hvplot were used to plot the Redlining overlay to see if there was any visual difference or benefits to using one type over the other, but it seems that both share similar conclusions. The benefit to hvplot in particular is that a mouse could be hovered over the plot and it will show coordinates, so this helpful to describe where exactly a particular feature is. Matplotlib and hvplot also have different parameters that can be used, so the use of a certain type of module for plotting is personal preference and if a certain parameter was essential. For this case study, both were done out of curiousity being new to the Earth Data Analysis field. While this (hvplot) and the matplotlib plot could be visually be compared to the actual redlining grade by looking back at the Redlining Map, it would be less prone to human error and less time consuming by caluculating zoning statistics and fitting a model to the data. This will share more specific information as to the relationship between redlining grade/category and green space (via NDVI).

# Calculate and Plot Zoning Stats Part 1 of 3

# Define new variable for denver redlining mask, using regionmask
philadelphia_redlining_mask = regionmask.mask_geopandas(
    # Put gdf in same CRS as raster
    philadelphia_redlining_gdf.to_crs(philadelphia_ndvi_da.rio.crs),
    # x and y coordinates from raster data 
    philadelphia_ndvi_da.x, philadelphia_ndvi_da.y,
    # The regions do not overlap
    overlap=False,
    # We're not using geographic coordinates
    wrap_lon=False,
)

# Plot the redlining mask to confirm regionmask worked
philadelphia_redlining_mask.plot(
    cbar_kwargs={"label": "Region"},
    )
plt.gca().set(
    title='Philadelphia Redlining Mask', 
    xlabel='',
    ylabel='',
    xticks=[], 
    yticks=[]
)
plt.show()

# Calculate and Plot Zoning Stats Part 2 of 3

# Calculate NDVI stats for each redlining zone
philadelphia_ndvi_stats = zonal_stats(
    philadelphia_redlining_mask,
    philadelphia_ndvi_da
)

# Calculate and Plot Zoning Stats Part 3 of 3
# Plot the regional statistics

# Merge the NDVI stats with redlining geometry into one `GeoDataFrame`
philadelphia_redlining_ndvi_gdf = philadelphia_redlining_gdf.merge(
    philadelphia_ndvi_stats.set_index('zone'),
    left_index=True, right_index=True
)
# Change grade to ordered Categorical for plotting
philadelphia_redlining_ndvi_gdf.grade = pd.Categorical(
    philadelphia_redlining_ndvi_gdf.grade,
    ordered=True,
    categories=['A', 'B', 'C', 'D']
)

# Drop rows with NA grades
philadelphia_redlining_ndvi_gdf = philadelphia_redlining_ndvi_gdf.dropna()

# Plot NDVI and redlining grade in linked subplots
(
    philadelphia_redlining_ndvi_gdf.hvplot(
        c='mean',
        geo=True,
        xaxis='Longitude', 
        yaxis='Latitude',
        title='Philadelphia Mean NDVI',
        clabel='Mean NDVI',
        cmap='Greens'
    )
    +
    philadelphia_redlining_ndvi_gdf.hvplot(
        c='grade',
        geo=True,
        xaxis='Longitude', 
        yaxis='Latitude',
        title='Philadelphia Redlining Grades',
        cmap='cet_diverging_bwr_20_95_c54'

    )

)

Low mean NDVI appears to visually be where the lowest or worst ranked redlining grade/ category is - there is a connection bewteen vegetation health (urban green spaces) and redlining grade/category

The lowest mean NDVI appears white or extremely light green on the plot and is located in the mostly in the south and east compared to the mid to high mean NDVI which is green to dark green and is located in the north and west. The plot for the Philadelphia redlining grades has a similar pattern in regards to the locations of the two worse/worst grades and the better/best grades. On this redlining grade plot the C (pink) and D(dark red) grades are in the south and east and the A (dark purple) and B (light purple) in the north and west. This similarity in location of the low mean NDVI also being where the C and D redlining grades are suggests a relationship between them - that the 2 lower redlining grades have a relationship to the lack of urban green spaces in the 21st century.

The location itself of where the lower 2 redlining grades are and where the low mean NDVI is, are in the oldest and most central part of the city which are in the south and east. This is the 'downtown' of Philadelphia and also where the 7th ward was when W.E.B. du Bois wrote about it. While the population of this ward was later displaced, the ward and its history serve as an example of downtown Philadelphia and the effects on this population. These two plots side by side visually tell a story or reinforce what was learned in the 'Site Description' that redlining has made a lasting impact on urban green spaces or the lack there of in these lower graded redlinined areas in the 21st century.

# Fit a model Part 1 of 3
# Fit a model

# Convert categories to numbers
philadelphia_redlining_ndvi_gdf['grade_codes'] = (
    philadelphia_redlining_ndvi_gdf.grade.cat.codes
)

# Fit model
philadelphia_tree_classifier = DecisionTreeClassifier(max_depth=2).fit(
    philadelphia_redlining_ndvi_gdf[['mean']],
    philadelphia_redlining_ndvi_gdf[['grade_codes']],
)

# Visualize tree
plot_tree(philadelphia_tree_classifier)
plt.show()

# Fit a model Part 2 of 3
# Plot the model results

# Predict grades for each region
philadelphia_redlining_ndvi_gdf ['predictions'] = (
    philadelphia_tree_classifier.predict(
        philadelphia_redlining_ndvi_gdf[['mean']])
)

# Subtract actual grades from predicted grades
philadelphia_redlining_ndvi_gdf['error'] = (
    philadelphia_redlining_ndvi_gdf ['predictions']
    - philadelphia_redlining_ndvi_gdf ['grade_codes']
)

# Plot the calculated prediction errors as a chloropleth
philadelphia_redlining_ndvi_gdf.hvplot(
    c='error',
    geo=True,
    xaxis='Longitude', 
    yaxis='Latitude',
    clabel='Predicted Grades Error',
    title='Philadelphia Calculated Prediction Errors',
    )

# Fit a model Part 3 of 3

# Evaluate the model with cross-validation
cross_val_score(
    DecisionTreeClassifier(max_depth=2),
    philadelphia_redlining_ndvi_gdf[['mean']],
    philadelphia_redlining_ndvi_gdf[['grade_codes']],
    cv=3
)

array([0.21428571, 0.2962963 , 0.25925926])

Model Description and Evaluation

For this Model, the 'scikit_learn' package was employed to use the class 'DecisionTreeClassifeir' to fit a Tree Model to the previously made 'philadelphia_redlining_ndvi_gdf' variable columns 'mean' (NDVI) and 'grade_codes' (converted redlining categories to numbers). After, this model was plotted using the function 'plot_tree' and module 'matplotlib.pyplot', revealed a simple model with only 4 recatangles (equal to the number of redlining categories) which is appropriate. This model with only 4 recatnagles/ thresholds was set by setting the parameter of 'max_depth' to 2, which limits the depth and overfitting of the model. This is an advantage of using this model. It is difficult to work with other possible choices because this particular dataset is so small that it's hard to see the diminishing point of return with other models.

Model results were further plotted to predict (Redlining) grades for each zonal region using the '.predict()' method of the 'DecisionTreeClassifier' class. By using 'hvplot' to plot the error of the predictions shows that there are very few (5) regions where the redlining grade/category was either off by positive 1 or 2 in that the redlining grade/category is actualy 1 or 2 better than predicted (i.e. a category B (grade_code=1) was predicted as C (grade_code=2) -
off by positive 1, or a category B was predicted as D - off by positive 2). However, there were many regions that were off by -1 or -2 in that the redlining grade/category is actualy 1 or 2 worse than predicted (i.e. a category C (grade_code=2) was predicted as B (grade_code=1) off by negative 1, or a category C (grade_code=2) was predicted as A (grade_code=0) off by negative 2). (The redlining map key and description is at the end of this section to better understand the redlining categories.)

Next, the model was evaluated by using the function 'cross_val_scores' to validate the accuracy of the model. This cross validation is beneficial because it also is good at identifying overfitting; however, because this datset is so small it can result in splits that are too small to remain meaningful, or not have all the categories. The evaluation showed values of [0.21428571, 0.2962963 , 0.25925926] which some are percentages greater than 25%, this means there is a relationship between mean NDVI and the redlining category, but one cross valuation value is below 25% and would not be better than random.

While the evaluation of the model doesn't present values that would indicate a strong relationship between urban grenspaces and redlining grades, there is a relationship, and the plot of the calculated predicted errors is visually telling. The predicted grades error is off by -2 in areas that are mostly redlining grade/category 'Hazordous' (D)(grade_code=3) and 'Definitely Declining' (C)(grade_code=2), meaning that the areas were predicted to be a better grade/category than they are in actuality. Because most of the calculated prediction errors are of this kind (either -1 or -2 - actual grade is worse than predicted) means that vegetation health in neighborhoods is related to its redlining grade/category in that there is better vegetation health, or more ubran green spaces, than would be expected. This could be for a few reasons, there are many urban greespaces and parks that are central to Philadelphia that are in those areas that are redlining grades/categories D or C which has to do with how the city was first built with the 5 public squares and then later during the City Beautiful Movement and the New Deal Era the expansion of urban green spaces for parkways and and parks around bodies of water (which almost completely surrounds the city). Also, the oldest and more dense part of the city is so built out or developed that at some point it becomes dificult to creat urban green spaces when all of the land is already used.

So, while the model evaluation doesn't have very high percentages, the plot itself does show that there is a relationship between vegetation health and redlining grade/ category.

Redlining map of Philadelphia - Key, PA courtesy of Mapping Inequality ( Nelson and Winling (2023) )

Conclusion

There were many images and plots created throughout this portfolio post to show the possible relationship between redlining grades and urban green spaces.

• Redlining Outline on Esri Imagery Tiles - oriented the viewer to current urban green spaces in relation to the redlining outline of Philadelphia.

• True Color and CIR Images - While the True Color Image visually showed the larger water ways and urban green spaces, the CIR Image made the smaller green spaces and water ways easier to see due to the contrast. The larger urban green spaces are near bodies of water near the edges of the city compared to the smaller green spaces which are more central to downtown and dispersed.

• NDVI with Redlining Overlays (matplotlib and hvplot) - visually could see a difference in green spaces and how the larger green spaces are like physical boundaries between better and worse graded Redlining areas. However, this analysis is prone to human error and this relationship could be further analyzed by other methods.

• Zonal Stats Side by side Mean NDVI and Redlining Grade - this side by side was helpful to see the relationship between low Mean NDVI and the worse 2 redlining grades meaning that redlining has made a lasting impact on urban green spaces or the lack thereof in these lower graded redlined areas in the 21st century.

• Calculated Prediction Errors and Model Evaluation - the model evaluation didn't have very high percentages, but most were better than random. Also, the plot itself did show that there is a relationship between vegetation health and redlining grade/ category - the mean NDVI would've predicted better redlining grades in the south and east where the redlining grades are actually the worst grades.

These plots and images built upon one another to show that there is a relationship between redlining grade and urban green spaces in that redlining grades C and D have lower mean NDVI as well as a visual lack of urban green spaces. Because urban green spaces were a large push of the City Beautiful Movement which was about increasing the city's prestige, the green spaces and civic institutions that were created involved slum clearing and displacing current residents to 'beautify' the city at expense of people, and communities, particularly those of color.

One can even see how Evenezer Howard's 'Garden City' description - "central park containing important public buildings and surrounded by a 'Crystal Palace' ring of retail stores ... encircled by a permanent agricultural greenbelt of some 5,000 acres" (Howard, 328), can be seen in Philadelphia with its Benjamin Franklin Parkway, and that area being lined with important civic and public buildings that are seen as beautiful and adding prestige and grandeur. What should be noted that was not included in the description of the 'Garden City' is how everyone, all people have history and have culture and roots in the current built and physical environment. In Philadelphia, communities and people of color, particularly Black communities were to be excluded from this idea beautifying of the city entirely and displaced to a further away area so this centralized beauty and grandeur would not be tainted by them. This is a prime example of structural racism.

Racism and utopian ideals that had force through the mechanism of redlining have influenced the distribution of urban green space in Philadelphia in the 21st century. While there has been work done to undo these redlining policies and create hopefully someday equal access and use to urban green space with the creation of new green spaces in these C and D redlined areas, its impact is still with us today.