Unlocking Spatial Insights: Loading County Parcel Data from Esri Feature Services into Geopandas for Spatial Analysis

Unlock the power of real estate analysis with this series! Learn to connect to ESRI Feature Services, load county parcel data into Geopandas, and explore insights without an ArcGIS License. Master these skills for enhanced expertise and discover hidden opportunities in upcoming articles.

Getting the Data:

In this article lets assume we are going to grab Parcel Data for Bastrop County Texas through their interactive Address Lookup Map:

Spotting the Land Parcel layer in the list, the next step is accessing it. Simply click on the three dots and navigate to ‘Show Item Details.

If you click on this link, it will redirect you to a specific layer item within the Bastrop County GIS ArcGIS Portal. On this page, you’ll find a wealth of information about the data, such as Terms of Use, Metadata, Updated Date times, Descriptions, and a Data Table Preview, among other details.

For our specific needs, focus on locating the URL, which is typically positioned towards the bottom right of the page.

The URL you find leads to the ESRI Feature Server, where the data for this specific layer is stored and can be accessed. This URL allows us to engage in various activities, such as accessing individual layers within the FeatureServer (though, in our instance, there is only one layer). We can perform tasks like exports and filtering. We can do more with ArcGIS GeoProcessing Credits but it’s worth noting that, for our current use case, no licensing is required.

Upon following the URL, we will arrive at the ParcelARI Feature Server, specifically within the ArcGIS Rest Services Directory. Although the layout may not be as visually appealing as previous pages, this section is crucial. Here, we directly access the data pertinent to our needs, specifically the bastrop_publish.bastrop_gis_sys.ParcelARI (0) layer.

Upon reaching this page, we are now examining the Feature Layer for our Parcels, which will be the focus of our queries. To simplify, let’s highlight the key elements on this page. We have the Layer Name, crucial for identification. The Type is noteworthy, particularly for our goal of accessing individual features, as it indicates that it is a Feature Layer. The geometry type is Polygon, signifying the shape of the parcels.

For practical considerations, take note that the maximum number of features (parcels) that can be queried in a single request is 1000 (This is important to remember). If more are needed, requests to their GIS can be made to adjust this limit, subject to their configuration.

The Geospatial Extent provides a view of the Earth’s surface covered by this data. The coordinates may seem unfamiliar but are based on theSpatial Reference 102739 (2277), a specific system for the Bastrop County Region. If you desire further information on Spatial References and Coordinate systems, feel free to ask.

Beyond this, the page includes generic details, and there is information about the fields within our data, offering a preview of the available information.

For right now this is all we care about — go ahead and copy the URL for this page — this is our Feature Layer and this is how we will programmatically access the data.

Using Python to Interact with the data:

So before anything in your Python environment (We are using Python 3.10) pip install geopandas and pip install requests.

These installations will grant you access to Geopandas, enabling enjoyable spatial processing and analysis, and Requests, a powerful library for making HTTP requests to URLs. In our case, we’ll leverage Requests to interact with ESRI REST Services for our Feature Layer. Now, let’s get started with the rest of the script

1. Requesting Features

import requests

# URL to the Bastrop County Parcel Feature Layer endpoint
feature_layer_url = "https://services3.arcgis.com/wdTkTU0MdZbNBEZy/ArcGIS/rest/services/ParcelARI/FeatureServer/0"

# Generic Parameters for right now - nothing fancy
params = {
    'where': '1=1',  # Query to retrieve all features
    'outFields': '*',  # * gets all fields
    'returnGeometry': True,  # We want Geometry to do spatial analysis
    'f': 'json'  # Response format as JSON
}

# Make the request
response = requests.get(feature_layer_url + '/query', params=params)

This makes the request against our Feature Layer and returns everything it possibly can (Which in this case is 1000 features)

2. Processing the Data into a Geodataframe

Okay so now we have our response we want to ensure it was successful and then we want to access the feature geometries and attributes and set our polygon geometries and create fields for our attributes:

import geopandas as gpd

# Check if the request was successful 
if response.status_code == 200:
    # Parse the JSON response
    data = response.json()

    # Access the features
    features = data.get("features", [])

    # Convert Esri JSON to GeoJSON for each feature
    for feature in features:
        geometry_esri = feature.get("geometry")
        if geometry_esri:
            geometry_geojson = {
                "type": "Polygon",
                "coordinates": geometry_esri.get("rings")
            }
            feature['geometry'] = geometry_geojson
            feature["properties"] = feature.get("attributes", {})

    # Create our GeoDataFrame
    parcel_gdf = gpd.GeoDataFrame.from_features(features, crs="epsg:2277")

    # Display the GeoDataFrame
    print(parcel_gdf.head())

else:
    print(f"Error: {response.status_code} - {response.text}")

Geopandas simplifies much of the data processing for us. Our primary task involves setting the correct geometry feature[‘geometry’] and property feature[‘properties’] to ensure the proper structuring of our data as a GeoJSON, ready to be loaded into a Geopandas dataframe. Furthermore, we establish our coordinate reference system (CRS) as EPSG:2277, based on the earlier observation in the metadata indicating that our data is in Spatial Reference 102739 (2277).

If we export this data to a shapefile with something like

# Export the GeoDataFrame to a Shapefile
parcel_gdf.to_file("parcel.shp", driver="ESRI Shapefile")

We can view it in an opensource GIS Software like QGIS or visualize with other Python packages — we can see the data we have exported:

Okay so this is very basic and we are still limited to 1000 records but hey its a start right?

3. Filtering Fun

We can filter and manipulate the data in various ways, but for simplicity, let’s focus on 2 criteria:

Filter by Land Value > $700,000

# Specify your query parameters with an additional filter for land_val
params = {
    'where': 'land_val > 700000',  # Filter for parcels with land_val greater than 700,000
    'outFields': '*',  # * gets all fields
    'returnGeometry': True,  # We want Geometry to do spatial analysis
    'f': 'json'  # Response format as JSON
}

# Make the request
response = requests.get(feature_layer_url + '/query', params=params)

It’s as simple as changing our where parameter to include a simple greater than expression

Using the previous code we made to create the geodataframe and output the geodataframe to shapefile we can see our $700,000 + parcels

Lat Lon Radius Search

from pyproj import Proj, transform

# Specify the specific latitude and longitude
lat, lon = 30.119838750410207, -97.31691735686985 

# Define the target CRS (EPSG:2277)
target_crs = Proj(init='epsg:2277')

# Transform the input coordinates to the target CRS
lon, lat = transform(Proj(init='epsg:4326'), target_crs, lon, lat)

# Specify your query parameters with filters for a specific location and 1-mile radius
params = {
    'geometry': f'{lon},{lat}',  # Set the center point
    'geometryType': 'esriGeometryPoint',
    'spatialRel': 'esriSpatialRelIntersects',
    'distance': 1,  # Set the radius in miles
    'units': 'esriSRUnit_StatuteMile',
    'outFields': '*',  # * gets all fields
    'returnGeometry': True,  # We want Geometry to do spatial analysis
    'f': 'json',  # Response format as JSON
    'outSR': 2277  # Specify the target CRS (EPSG:2277)
}

# Make the request
response = requests.get(feature_layer_url + '/query', params=params)

Take a look at this interesting capability. We want to use our input Latitude and Longitude to filter Parcels within a 1 mile radius. Given that our Parcel data uses CRS 2277 while our input coordinates are in the more common 4326 format, we use Pyproj to transform our input coordinates. Following that, we specify that the query should execute an Intersection, capturing everything within our defined 1-mile radius. It's a straightforward process that efficiently gets the job done.

Using the previous code we made to create the geodataframe and output the geodataframe to shapefile we can see our Radius results

To conclude

That provides a brief overview of ESRI Rest capabilities, showcasing that interaction with the data doesn’t necessarily require an ArcGIS License. If you’re considering exporting an entire county’s worth of data, it’s advisable to check for export options on the County GIS Site or explore potential availability through a request contact form. This glimpse offers a simple approach to engaging with the GIS data available for your county. Stay tuned for upcoming articles as I plan to delve deeper into the subject of county GIS data