Seniors at Risk: Using Spatial Analysis to Identify Pharmacy Deserts

Seniors are especially vulnerable during the COVID-19 pandemic.

At makepath, we build and use open source spatial analysis tools such as Xarray-Spatial and Datashader.

In this post, we show you how to use those tools to identify and quantify pharmacy deserts where seniors are at especially high risk.

We will walk through our analysis to show who is at risk. For readers who want to dive deeper, we provide a link to an example notebook showing the implementation using Xarray-Spatial and Datashader.

With analysis results in hand, we called people who live in these areas to get their perspectives and put faces and voices to the problem of pharmacy deserts.

You can apply these free and open source spatial analysis tools to your own problem domains while using this case study and accompanying notebook as a starting place.

At makepath we specialize in custom implementations to tackle big data questions, especially those with a focus on geospatial analysis. If you have any questions, please email us at

What is a Pharmacy Desert?

A pharmacy desert is an area with reduced access to prescription and over-the-counter medications. Older residents are more vulnerable to pharmacy deserts because they are more dependent on these essential supplies.

Understanding pharmacy deserts is important to businesses looking to expand and to government agencies looking to guide development, provide services, and save lives.

Using spatial analytics and publicly available data, we identified 10 counties with high percentages of residents 65 and older who live in pharmacy deserts.

These counties include:

  • Hooker County, Nebraska
  • Catron County, New Mexico
  • Harding County, New Mexico
  • Wheeler County, Oregon
  • Jeff Davis County, Texas
  • Prairie County, Montana
  • Edwards County, Texas
  • Real County, Texas
  • Idaho County, Idaho
  • Valley County, Idaho

We identified these counties using the raster-based spatial analysis tools within Xarray-Spatial. Our conclusion is that in counties like Hooker, Catron and Harding, pharmacies have an opportunity to serve these communities.

In counties like Hooker, Catron and Harding, pharmacies have an opportunity to serve these communities.

Let’s now understand how we arrived at these 10 counties. The following analysis can be performed on consumer-grade computers using free and open source tools with publicly available data.

Step 0. Get Publicly Available Data

Our analysis of pharmacy deserts starts by obtaining pharmacy locations from the Department of Homeland Security (DHS). The DHS updates pharmacy locations regularly with precise latitude/longitude locations via their ArcGIS Online organization.

With pharmacy locations downloaded, we now need age information. ArcGIS Online also makes census information available pre-joined with age and sex demographics. Here we are able to obtain census blockgroups with our desired attributes as a shapefile. Let’s also grab county boundaries from the US Census Bureau.

In our analysis notebook, we load data into memory using pandas (pharmacy locations), and geopandas (block groups, counties).

pharmacies_df = pandas.read_csv("pharmacy_locations.csv"))
blockgroup_df = GeoDataFrame(geopandas.read_file("bg_wgs84.shp"))
counties_df = GeoDataFrame(geopandas.read_file("counties_wgs84.shp"))

These three layers (pharmacies, blockgroup ages, and counties) are all we need to begin an exploratory spatial analysis of pharmacy deserts. We will be performing raster-based analysis to align data to a regular grid overlay with layers on top of one another.

Step 1. Define Study Area

Map in Datashader
1600×800 pixel mask of lower-48 United States Study Area

We start our analysis by defining the study area. We’ll spatially align our data to this study area’s resolution. The precision of our analysis is controlled by the resolution of this layer.

Step 2. Create Distance to Pharmacy Layer

Spatial Analysis map made with Datashader and Xarray-Spatial
Proximity grid to nearest pharmacy classified into 5 categories using natural breaks. Lighter color means farther away from pharmacy.
distance = proximity(pharmacy_raster, distance_metric='GREAT_CIRCLE')
classified_distance = natural_breaks(distance, k=5)

Using pharmacy latitude / longitude locations, we can generate a proximity grid which shows the distance of every pixel to the nearest pharmacy location. Since we are dealing with geographic coordinates, we use great circle distance to compute distance over a spheroid. Then we classify the distances into 5 categories using natural breaks.

Step 3. Create Percent Older than 65 Layer

Spatial Analysis map made with Datashader and Xarray-Spatial
Percent of residents over 65-years-old classified into 5 categories using natural breaks. Lighter means “great percentage of +65 residents”.
classified_age = natural_breaks(percent_pop_65_raster, k=5)

Using US Census Block Groups, we can create a raster based on the percentage of 65+ residents in the block group. We then take this age raster and classify it into 5 categories using natural breaks.

Step 4. Combine Layers to Find Areas of Interest

Spatial Analysis map made with Datashader and Xarray-Spatial
This binary raster (values either 0 or 1) shows areas of interest defined as locations of high percentages of senior citizen AND far from nearby pharmacy.
is_pharmacy_desert = binary(classified_distance, [5]) # only category 5
is_high_percent_above_65 = binary(classified_age, [5]) # only category 5
area_of_interest = is_pharmacy_desert * is_high_percent_above_65

With distance and age layers defined, we can create binary (0 / 1) layers for both and multiply to get our areas of interest. When we use the Xarray-spatial binary classification tool, we get back a new layer of only ones and zeros. When we multiply these two 0/1 layers together, all locations with a zero drop out and only locations with 1 in both age and distance categories remain. We’ve now been able to isolate our counties of interest: those that are far from pharmacies and have a high percentage of senior citizens.

Step 5. Summarize Areas of Interest by County

Spatial Analysis map made with Datashader and Xarray-Spatial
County boundaries from US Census Bureau which are used to summarize our pharmacy desert layer.
results = zonal_stats(counties_raster, area_of_interest, zonal_funcs)
results.nlargest(10, 'pharmacy_desert_score')

Zonal statistics allow us to summarize our Areas of Interest by county. The results are individual pharmacy_desert_scores for each county. We can then take the top 10 score which is our final list of counties to target our outreach efforts.

Jeff Davis0.47
Output table from Zonal Statistics operation which summarizes our pharmacy desert score raster for each county in the lower-48 United States. We then take the top 10 counties by pharmacy_desert_score.

Picking Up the Phone

We used tools available in Xarray-Spatial and Datashader to visualize and quantify pharmacy deserts by county, with a focus on areas comprised of older residents at risk from reduced access to medications.

Then it was time to pick up the phone. From this simple analysis, we set up phone calls with residents of these communities to learn more about their lives. These interviews helped us ground-truth the analysis and verify our results.

  • Jennifer Peters lives in Hooker County, Nebraska. She confirmed the challenging reality of living 72 miles from the nearest pharmacy. “Hooker County residents rely on delivery services for prescriptions, and connect on social media to share meds,” said Peters.
  • Scott Johnston is the County Assessor in Catron County, New Mexico. He shared that a large percentage of county residents are above 65 years-old and are negatively impacted by the two-hour round trip drive to get prescriptions filled.
  • Virginia Smith is the Senior Citizen Program Coordinator in Harding County, New Mexico. She described the stress of the three-hour round trip to fill prescriptions. Many Harding County residents do not have a credit card, which is necessary to get medications by mail.

This analysis highlights areas with higher percentages of senior citizen that exist in pharmacy deserts. Our analysis is limited to the contiguous United States, but the techniques can be easily applied to different study areas and topics.

You can find the full analysis in the Xarray-Spatial github repo.


Spatial analysis can be a powerful tool to answer real world questions. Using open source tools, we were able to quickly identify pharmacy deserts. We then confirmed our results by calling people and listening to their challenges.

This was an exploratory analysis. If we were to go further, our next step would be to do a more localized study of the problem areas we initially identified.

This would include zoomed-in maps and a more thorough site selection analysis for future pharmacy builds. We would enrich the layers with additional characteristics of the population and average distance to pharmacies overlaid with available commercial real estate parcels.

We would also conduct a larger number of interviews to ground-truth the conclusions with a larger degree of confidence.

Have you used open source tools to answer this kind of question? What has been your experience with pharmacy deserts or other similar problems? Do you need help developing custom GIS solutions to answer complex geospatial questions?

Please let us know in the comments, or reach out at

2 thoughts on “Seniors at Risk: Using Spatial Analysis to Identify Pharmacy Deserts

  1. Pingback: Spatial Conversations – MakePath & xarray-spatial – VerySpatial

  2. Pingback: In-Depth Video: Pharmacy Desert Identification with Open Source Spatial Analysis Tools - makepath

Leave a Reply

Your email address will not be published.