Redlining Analysis Reports

Redlining 93 Denver v1

Denver Redlining

Introduction

Redlining is a discriminatory practice in which financial services are systematically withheld from neighborhoods with significant racial and ethnic minority populations. This practice has been most prominent in the United States, predominantly targeting African-American communities. Common examples of redlining include the denial of credit, insurance, healthcare services, and the creation of food deserts in minority neighborhoods.

Redlining in the U.S. emerged within a broader context of racial segregation and discrimination, shaped by economic theories on race and property values promoted by Richard T. Ely’s Institute for Research in Land Economics. The National Housing Act of 1934 and the establishment of the Federal Housing Administration (FHA) institutionalized this practice, under the leadership of Homer Hoyt, who developed the first mortgage underwriting criteria. This federal policy accelerated the decline of minority neighborhoods by restricting access to mortgage capital, exacerbating residential segregation and contributing to urban decay.

Redlining in Denver, which began in the 1930s, was a discriminatory practice that systematically excluded minority communities, particularly African Americans, from access to financial services and housing opportunities. The Federal Housing Administration (FHA) created maps that marked neighborhoods deemed "high-risk" in red, discouraging banks from issuing mortgage loans in these areas.

One notable example is the Five Points neighborhood, which, despite the restrictions imposed by redlining, emerged as a vibrant cultural hub known as the "Harlem of the West." From the 1920s to the 1960s, Five Points thrived with a rich jazz scene and an active African American business community.

Although redlining was officially banned in 1968, its effects persisted, contributing to residential segregation and a lack of investment in certain neighborhoods. In recent years, Denver has implemented programs to address these historical inequalities. For example, in 2022, the city launched an initiative offering financial assistance for down payments to individuals affected by redlining between 1938 and 2000, aiming to promote equity in homeownership.

Additionally, studies have shown that areas historically impacted by redlining in Denver continue to face significant challenges in terms of health equity and economic opportunities, reflecting the long-term consequences of these discriminatory practices.

In summary, redlining in Denver had a profound and lasting impact on the city’s development, perpetuating segregation and limiting opportunities for minority communities. Current efforts seek to acknowledge and remedy these historical injustices to build a more inclusive and equitable city.

Data Citation

• "Redlining: How the Five Points Neighborhood Formed Amid Racist Practices in the 1930s"
• "Redlined in Denver Between 1938 and 2000? The City Might Help You Put a Down Payment on a New Home"
• "Denver’s Public Health and Environment Equity Report (2020)"

These sources detail the historical impact of redlining in Denver, the specific case of Five Points, and current efforts to address the effects of this discriminatory practice.

Description of the Resources

A multispectral image will be imported, product is HLSL30, we will use the green band.

Multispectral image: A set of images taken in different bands of the electromagnetic spectrum (e.g., visible bands such as red, green, blue and others in the near infrared). The revisit days is 8 of the sensor. This product provides images with a spatial resolution of 30 meters, adjusted to correct for bidirectional reflectance distribution function (BRDF) effects, allowing global observations of the Earth's surface. We will use the reflectance of the product.

HLSL30” product: Is a identifier of a dataset, sensor, processed image, or geospatial product generated by remote sensing or satellite image analysis.

Green band: The spectral channel corresponding to wavelengths associated with the color green, often used in remote sensing to assess vegetation, photosynthesis, and crop health. Spectral channel (Green): 0.53 – 0.59 µm.

Import Libraries

%store -r denver_redlining_gdf data_dir

import re # Use regular expressions to extract metadata

import earthaccess # Access NASA data from the cloud
import matplotlib.pyplot as plt  #Overlay raster and vector data
import numpy as np # Process bit-wise cloud mask
import pandas as pd # Group and aggregate
import rioxarray as rxr # Work with raster data
from rioxarray.merge import merge_arrays # Merge rasters

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

denver_results = earthaccess.search_data(
    short_name="HLSL30",
    bounding_box=tuple(denver_redlining_gdf.total_bounds),
    temporal=("2023-07-12"),
    count=1
)

denver_files = earthaccess.open(denver_results)

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

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

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

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
    """
    #Connet 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

process_image(denver_files[8], denver_redlining_gdf).plot()

<matplotlib.collections.QuadMesh at 0x20447652dd0>

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 significat bit first to match the user guide
            bitorder='little',
            # Expand the array in a new dimension
            axis=-1)
    )


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

    return cloud_mask

blue_da = process_image(denver_files[1], denver_redlining_gdf)
denver_cloud_mask = process_cloud_mask(
    denver_files[-1],
    denver_redlining_gdf,
    [1, 2, 3, 5])
#blue_da.where(denver_cloud_mask).plot()
    
  

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

denver_results = earthaccess.search_data(
    short_name="HLSL30",
    bounding_box=tuple(denver_redlining_gdf.total_bounds),
    temporal=("2023-07-12"),
)

# Open earthaccess results
denver_files = earthaccess.open(denver_results)

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

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

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

denver_files[8].full_name

'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T13SDD.2023193T173653.v2.0/HLS.L30.T13SDD.2023193T173653.v2.0.B03.tif'

# Compile a regular expression to search for metadata
uri_re = 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
uri_groups = [
    uri_re.search(denver_file.full_name).groupdict()
    for denver_file in denver_files]

# Create a DataFrame with the metadata
raster_df = pd.DataFrame(uri_groups)

# Add the File-like URI to the DataFrame
raster_df['file'] = denver_files

# Check the results
#raster_df

# Labels for each band to process
bands = {
    'B02': 'red',
    'B03': 'green',
    'B04': 'blue',
    'B05': 'nir'
}
# Initialize structure for saving images
denver_das = {band_name: [] for band_name in bands.values()}
for tile_id, tile_df in raster_df.groupby('tile_id'):
     # Load the cloud mask
    fmask_file = tile_df[tile_df.band_id=='Fmask'].file.values[0]
    cloud_mask = process_cloud_mask(
        fmask_file,
        denver_redlining_gdf,
        [1, 2, 3, 5])
    

    for band_id, row in tile_df.groupby('band_id'):
        if band_id in bands:
           band_name = bands[band_id]
           print(band_id, band_name)
           # Process band
           band_da = process_image(row.file.values[0], denver_redlining_gdf)

           # Mask band
           band_masked_da = band_da.where(cloud_mask)
           

           # Store the resulting DataArray ofr later
           denver_das[band_name].append(band_masked_da)
    break

B02 red
B03 green
B04 blue
B05 nir

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

# Plot a merged raster band
denver_merged_das['green'].plot(cmap='Greens', robust=True)

<matplotlib.collections.QuadMesh at 0x2044e2d6750>

Redlining 41 Zonal Stats

STEP 6: Calculate zonal statistics

In order to evaluate the connection between vegetation health and redlining, we need to summarize NDVI across the same geographic areas as we have redlining information.

First, import variables from previous notebooks:

store -r denver_redlining_gdf denver_ndvi_da

Try It: Import packages

Some packages are included that will help you calculate statistics for areas imported below. Add packages for:

1. Interactive plotting of tabular and vector data
2. Working with categorical data in DataFrames

import hvplot.pandas # Interactive plots with pandas
import pandas as pd # Ordered categorical data
import regionmask # Convert shapefile to mask
from xrspatial import zonal_stats # Calculate zonal statistics

Try It: Convert vector to raster

You can convert your vector data to a raster mask using the regionmask package. You will need to give regionmask the geographic coordinates of the grid you are using for this to work:

1. Replace gdf with your redlining GeoDataFrame.
2. Add code to put your GeoDataFrame in the same CRS as your raster data.
3. Replace x_coord and y_coord with the x and y coordinates from your raster data.

denver_ndvi_da.plot(robust=True)

<matplotlib.collections.QuadMesh at 0x1fa0800d590>

denver_redlining_gdf.plot()

<Axes: >

denver_redlining_mask = regionmask.mask_geopandas(
    denver_redlining_gdf.to_crs(denver_ndvi_da.rio.crs),
    denver_ndvi_da.x, denver_ndvi_da.y,
    # The regions do not overlap
    overlap=False,
    # We're not using geographic coordinates
    wrap_lon=False
)
denver_redlining_mask.plot()

<matplotlib.collections.QuadMesh at 0x1fa0a27b710>

Try It: Calculate zonal statistics

Calculate zonal status using the zonal_stats() function. To figure out which arguments it needs, use either the help() function in Python, or search the internet.

# help(zonal_stats)

# Calculate NDVI stats for each redlining zone
denver_ndvi_stats = zonal_stats(denver_redlining_mask, denver_ndvi_da)
# denver_ndvi_stats

Try It: Plot regional statistics

Plot the regional statistics:

1. Merge the NDVI values into the redlining GeoDataFrame.
2. Use the code template below to convert the grade column (str or object type) to an ordered pd.Categorical type. This will let you use ordered color maps with the grade data!
3. Drop all NA grade values.
4. Plot the NDVI and the redlining grade next to each other in linked subplots.

denver_redlining_gdf.head()

	area_id	city	state	city_surve	category	grade	label	residentia	commercial	industrial	fill	geometry
1103	6525	Denver	CO	True	Best	A	A1	True	False	False	#76a865	POLYGON ((-104.90814 39.74754, -104.90824 39.7...
1104	6529	Denver	CO	True	Best	A	A2	True	False	False	#76a865	POLYGON ((-104.92221 39.73851, -104.92228 39.7...
1105	6537	Denver	CO	True	Best	A	A3	True	False	False	#76a865	POLYGON ((-104.91281 39.72538, -104.91779 39.7...
1106	6536	Denver	CO	True	Best	A	A4	True	False	False	#76a865	POLYGON ((-104.96016 39.72913, -104.96016 39.7...
1107	6540	Denver	CO	True	Best	A	A5	True	False	False	#76a865	POLYGON ((-104.95956 39.72366, -104.95947 39.7...

# Merge the NDVI stats with redlining geometry into one `GeoDataFrame`
denver_redlining_ndvi_gdf = denver_redlining_gdf.merge(
    denver_ndvi_stats.set_index('zone'),
    left_index=True, right_index=True)

# Change grade to ordered Categorical for plotting
denver_redlining_ndvi_gdf.grade = pd.Categorical(
    denver_redlining_ndvi_gdf.grade,
    ordered=True,
    categories=['A', 'B', 'C', 'D']
)

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

# Plot NDVI and redlining grade in linked subplots
(
denver_redlining_ndvi_gdf.hvplot(
    title='NDVI',
    c='mean',
    geo=True,
    cmap='Greens')
+
denver_redlining_ndvi_gdf.hvplot(
    title='Denver Redlining',
    c='grade',
    geo=True,
    cmap='cet_cwr')    
)

store denver_redlining_ndvi_gdf

Stored 'denver_redlining_ndvi_gdf' (GeoDataFrame)

Redlining 42 Tree Model

STEP 7: Fit a model

One way to determine if redlining is related to NDVI is to see if we can correctly predict the redlining grade from the mean NDVI value. With 4 categories, we’d expect to be right only about 25% of the time if we guessed the redlining grade at random. Any accuracy greater than 25% indicates that there is a relationship between vegetation health and redlining.

To start out, we’ll fit a Decision Tree Classifier to the data. A decision tree is good at splitting data up into squares by setting thresholds. That makes sense for us here, because we’re looking for thresholds in the mean NDVI that indicate a particular redlining grade.

First, import variables from previous notebooks:

store -r denver_redlining_ndvi_gdf

Try It: Import packages

The cell below imports some functions and classes from the scikit-learn package to help you fit and evaluate a decision tree model on your data. You may need some additional packages later one. Make sure to import them here.

import hvplot.pandas
import matplotlib.pyplot as plt
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.model_selection import train_test_split, cross_val_score

As with all models, it is possible to overfit our Decision Tree Classifier by splitting the data into too many disconnected rectangles. We could theoretically get 100% accuracy this way, but drawing a rectangle for each individual data point. There are many ways to try to avoid overfitting. In this case, we can limit the depth of the decision tree to 2. This means we’ll be drawing 4 rectangles, the same as the number of categories we have.

Alternative methods of limiting overfitting include:

• Splitting the data into test and train groups – the overfitted model is unlikely to fit data it hasn’t seen. In this case, we have relatively little data compared to the number of categories, and so it is hard to evaluate a test/train split.
• Pruning the decision tree to maximize accuracy while minimizing complexity. scikit-learn will do this for you automatically. You can also fit the model at a variety of depths, and look for diminishing accuracy returns.

Try It: Fit a tree model

Replace predictor_variables and observed_values with the values you want to use in your model.

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

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

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

#plot the model result
denver_redlining_ndvi_gdf['predictions'] = (
    tree_classifier.predict(denver_redlining_ndvi_gdf[['mean']])
)
denver_redlining_ndvi_gdf['error'] = (
    denver_redlining_ndvi_gdf['predictions']
    - denver_redlining_ndvi_gdf['grade_codes']
)
denver_redlining_ndvi_gdf.hvplot(c='error', geo=True)

c:\Users\Anamaria\miniconda3\envs\earth-analytics-python\Lib\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)

Try It: Plot model results

Create a plot of the results by:

1. Predict grades for each region using the .predict() method of your DecisionTreeClassifier.
2. Subtract the actual grades from the predicted grades
3. Plot the calculated prediction errors as a chloropleth.

One method of evaluating your model’s accuracy is by cross-validation. This involves selecting some of your data at random, fitting the model, and then testing the model on a different group. Cross-validation gives you a range of potential accuracies using a subset of your data. It also has a couple of advantages, including:

1. It’s good at identifying overfitting, because it tests on a different set of data than it trains on.
2. You can use cross-validation on any model, unlike statistics like $p$-values and $R^2$ that you may have used in the past.

A disadvantage of cross-validation is that with smaller datasets like this one, it is easy to end up with splits that are too small to be meaningful, or don’t have all the categories.

Remember – anything above 25% is better than random!

Try It: Evaluate the model

Use cross-validation with the cross_val_score to evaluate your model. Start out with the 'balanced_accuracy' scoring method, and 4 cross-validation groups.

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

array([0.66666667, 0.64285714, 0.57142857])

Looking for an Extra Challenge?: Fit and evaluate an alternative model

Try out some other models and/or hyperparameters (e.g. changing the max_depth). What do you notice?

print(denver_redlining_ndvi_gdf.columns)

Index(['area_id', 'city', 'state', 'city_surve', 'category', 'grade', 'label',
       'residentia', 'commercial', 'industrial', 'fill', 'geometry', 'mean',
       'max', 'min', 'sum', 'std', 'var', 'count', 'grade_codes',
       'predictions', 'error'],
      dtype='object')

# Try another model
from sklearn.linear_model import LogisticRegression

logistic_model = LogisticRegression(max_iter=1000).fit(
    denver_redlining_ndvi_gdf[['mean']],
    denver_redlining_ndvi_gdf.grade_codes,
)

accuracy = logistic_model.score(
    denver_redlining_ndvi_gdf[['mean']],
    denver_redlining_ndvi_gdf.grade_codes,
)
print(f"Accuracy del modelo de regresión logística: {accuracy}")

Accuracy del modelo de regresión logística: 0.3953488372093023

print(denver_redlining_ndvi_gdf['grade_codes'].value_counts())

grade_codes
2    17
1    10
3    10
0     6
Name: count, dtype: int64

Reflect and Respond

Practice writing about your model. In a few sentences, explain your methods, including some advantages and disadvantages of the choice. Then, report back on your results. Does your model indicate that vegetation health in a neighborhood is related to its redlining grade?

The Logistic Regression Model:

The Logistic Regression model is a solid choice for binary or multiclass classification problems, especially when the relationships between the independent variables and the probability of the classes are approximately linear. However, its suitability depends on the specific context of your data and objectives. Analysis of Class Distribution

Majority Class

Class 2 represents 17 out of a total of 43 instances: Majority class ratio =17/43=39.5%

Conclusión

The model is probably biased towards the majority class (2) and is not learning useful patterns to distinguish between the smaller classes (0, 1, 3).