2024 New Jersey UAP Flap Location Visualization Code (Python)

import geopandas as gpd
import matplotlib.pyplot as plt
from geopy.geocoders import Nominatim
from shapely.geometry import Point

# List of locations in New Jersey
locations = [
    "Hillsboro, NJ", "Morristown, NJ", "Morris County, NJ", "Union County, NJ", "Scotch Plains, NJ",
    "Westfield, NJ", "Cranford, NJ", "Garwood, NJ", "North Plainfield, NJ", "Rutherford, NJ",
    "Sparta, NJ", "Sussex County, NJ", "Florham Park, NJ", "Burlington County, NJ", "Somerset County, NJ",
    "Marlton, NJ", "Glassboro, NJ", "Flemington, Hunterdon County, NJ", "Branchburg, NJ", "Readington, NJ",
    "Monmouth County, NJ", "East Brunswick, NJ", "Passaic County, NJ", "Warren County, NJ", "Wyckoff, NJ",
    "Ocean County, NJ", "Old Bridge, NJ", "Hazlet, NJ", "Matawan, NJ"
]

# Geolocate using Nominatim
geolocator = Nominatim(user_agent="nj_location_mapper")
geo_data = []
for location in locations:
    try:
        loc = geolocator.geocode(location)
        if loc:
            geo_data.append((location, loc.latitude, loc.longitude))
    except:
        geo_data.append((location, None, None))

# Create a GeoDataFrame
geometry = [Point(xy) for _, lat, lon in geo_data if lat and lon for xy in [(lon, lat)]]
valid_geo_data = [row for row in geo_data if row[1] and row[2]]
gdf = gpd.GeoDataFrame(valid_geo_data, columns=["Location", "Latitude", "Longitude"], geometry=geometry, crs="EPSG:4326")

# Load a New Jersey county shapefile, taken from here: https://catalog.data.gov/dataset/tiger-line-shapefile-current-state-new-jersey-county-subdivision
nj_counties_path = "./tl_2023_34_cousub/tl_2023_34_cousub.shp"  # Replace with the path to your shapefile
nj_counties = gpd.read_file(nj_counties_path)

# Plot New Jersey counties and locations
fig, ax = plt.subplots(figsize=(12, 10))
nj_counties.boundary.plot(ax=ax, linewidth=0.5, color='black')
nj_counties.plot(ax=ax, color="lightgray", edgecolor="gray", alpha=0.05)
gdf.plot(ax=ax, color='red', markersize=50)

# Annotate locations
#for x, y, label in zip(gdf.geometry.x, gdf.geometry.y, gdf['Location']):
#    ax.text(x, y, label, fontsize=8, ha='right')

plt.title("New Jersey Counties with Locations", fontsize=14)
plt.show()