# Ignore SettingWithCopyWarning warnings
warnings.filterwarnings('ignore', category=pd.errors.SettingWithCopyWarning)
# -----------------------------------------------------------------------------
# 1. LOAD DATA
# -----------------------------------------------------------------------------
geojson_file = "limites-administratives-agglomeration-nad83.geojson"
csv_file = "criminal_acts.csv"
# Load administrative boundary
try:
geo_data = gpd.read_file(geojson_file)
except Exception as e:
print(f"Error loading GeoJSON file: {e}")
exit()
# Load crime data
try:
criminal_acts_df = pd.read_csv(csv_file)
except Exception as e:
print(f"Error loading CSV file: {e}")
exit()
# -----------------------------------------------------------------------------
# 2. FILTER CRIME DATA (Vol de vĆ©hicule Ć moteur)
# -----------------------------------------------------------------------------
# Keep only rows with valid lat/lon
criminal_acts_df = criminal_acts_df.dropna(subset=["LATITUDE", "LONGITUDE"])
# Filter for "Vol de vĆ©hicule Ć moteur"
car_thefts_df = criminal_acts_df.query("CATEGORIE == 'Vol de vĆ©hicule Ć moteur'")
# -----------------------------------------------------------------------------
# 3. CREATE GEODATAFRAMES & REPROJECT
# -----------------------------------------------------------------------------
# Convert boundary to WGS84, then to EPSG:32198 (Lambert projection for Quebec)
geo_data = geo_data.to_crs(epsg=4326).to_crs(epsg=32198)
# Convert car thefts to a GeoDataFrame in EPSG:4326, then reproject to EPSG:32198
car_thefts_gdf = gpd.GeoDataFrame(
car_thefts_df,
geometry=gpd.points_from_xy(car_thefts_df["LONGITUDE"], car_thefts_df["LATITUDE"]),
crs="EPSG:4326"
).to_crs(epsg=32198)
# -----------------------------------------------------------------------------
# 4. CREATE A 1 KM x 1 KM FISHNET (GRID)
# -----------------------------------------------------------------------------
# Determine bounding box for the entire administrative boundary
minx, miny, maxx, maxy = geo_data.total_bounds
# Grid size (1 km = 1000 meters in Lambert projection)
grid_size = 1000.0
# Generate a list of polygons (1x1 km) covering the bounding box
polygons = []
x = minx
while x < maxx:
y = miny
x_right = min(x + grid_size, maxx)
while y < maxy:
y_top = min(y + grid_size, maxy)
polygons.append(box(x, y, x_right, y_top))
y += grid_size
x += grid_size
grid_gdf = gpd.GeoDataFrame(geometry=polygons, crs=geo_data.crs)
# (Optional) Intersect grid with boundary if you only want squares inside the admin area
grid_gdf = gpd.overlay(grid_gdf, geo_data, how='intersection')
# -----------------------------------------------------------------------------
# 5. SPATIAL JOIN: COUNT THEFTS PER 1KM GRID CELL
# -----------------------------------------------------------------------------
# Each theft point will be matched to the cell it falls within
joined = gpd.sjoin(car_thefts_gdf, grid_gdf, how='left', predicate='within')
# Count how many thefts per cell (index_right is the cell's index in the grid)
counts_per_cell = joined.groupby("index_right").size().rename("theft_count")
# Attach that count to the grid
grid_gdf["theft_count"] = grid_gdf.index.map(counts_per_cell).fillna(0).astype(int)
# -----------------------------------------------------------------------------
# 6. IDENTIFY THE SQUARE WITH THE HIGHEST THEFT COUNT
# -----------------------------------------------------------------------------
max_cell_index = grid_gdf["theft_count"].idxmax()
max_theft_count = grid_gdf.loc[max_cell_index, "theft_count"]
max_cell = grid_gdf.loc[[max_cell_index]] # as a GeoDataFrame
# -----------------------------------------------------------------------------
# 7. OBTAIN LATITUDE/LONGITUDE OF THIS SQUARE'S CENTROID
# -----------------------------------------------------------------------------
# First compute the centroid in the projected CRS (EPSG:32198)
max_cell_proj = max_cell.to_crs(epsg=32198)
centroid_proj = max_cell_proj.geometry.centroid.iloc[0]
# Convert that centroid to WGS84
centroid_wgs84 = gpd.GeoSeries([centroid_proj], crs="EPSG:32198").to_crs(epsg=4326)
centroid_lat = centroid_wgs84.geometry.iloc[0].y
centroid_lon = centroid_wgs84.geometry.iloc[0].x
print(f"Highest theft count: {max_theft_count}")
print(f"Centroid of the 1 kmĀ² cell in WGS84: Lat={centroid_lat:.6f}, Lon={centroid_lon:.6f}")
# -----------------------------------------------------------------------------
# 8. PLOT THE RESULTS: SHOW GRID + THEFT POINTS + RED OUTLINE AROUND MAX CELL
# -----------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the administrative boundary (no label to avoid PatchCollection in legend)
geo_data.plot(ax=ax, color="white", edgecolor="black", alpha=0.5)
# Plot the grid, colored by theft_count
grid_plot = grid_gdf.plot(
ax=ax,
column="theft_count",
cmap="Reds",
legend=True, # We'll create a colorbar
alpha=0.6,
edgecolor="grey",
legend_kwds={"label": "Number of Vehicle Thefts per kmĀ²"}
)
# Plot the theft points (no label to avoid PatchCollection legend warning)
car_thefts_gdf.plot(
ax=ax,
color="blue",
markersize=5,
alpha=0.5
)
# Highlight the top cell with a red outline (no fill)
max_cell.plot(
ax=ax,
facecolor="none",
edgecolor="red",
linewidth=3
)
ax.set_title("1 kmĀ² Grid with the Highest 'Vol de vĆ©hicule Ć moteur' Incidents")
ax.set_xlabel("X (meters, EPSG:32198)")
ax.set_ylabel("Y (meters, EPSG:32198)")
# -----------------------------------------------------------------------------
# CREATE A MANUAL LEGEND FOR NON-COLORBAR ITEMS
# -----------------------------------------------------------------------------
# 1) Boundary - white fill, black edge
boundary_patch = mpatches.Patch(facecolor='white', edgecolor='black', label='Boundary')
# 2) Theft incidents - blue marker
theft_marker = mlines.Line2D([], [], color='blue', marker='o', linestyle='None',
markersize=5, label='Car Theft Incidents')
# 3) Max cell highlight - red outline
max_patch = mpatches.Patch(facecolor='none', edgecolor='red', linewidth=2,
label=f'Highest count (N={max_theft_count})')
ax.legend(handles=[boundary_patch, theft_marker, max_patch], loc="upper right")
plt.tight_layout()
plt.show()
