geopandas

# Read various formats
gdf = gpd.read_file("data.shp")       # Shapefile
gdf = gpd.read_file("data.geojson")   # GeoJSON
gdf = gpd.read_file("data.gpkg")      # GeoPackage
 
# Read with spatial filter (faster for large files)
gdf = gpd.read_file("data.gpkg", bbox=(xmin, ymin, xmax, ymax))
 
# Write to file
gdf.to_file("output.gpkg")
gdf.to_file("output.geojson", driver="GeoJSON")
 
# PostGIS database
from sqlalchemy import create_engine
engine = create_engine("postgresql://user:pass@localhost/db")
gdf = gpd.read_postgis("SELECT * FROM table", con=engine, geom_col='geom')

# Check CRS
print(gdf.crs)
 
# Set CRS (when metadata missing)
gdf = gdf.set_crs("EPSG:4326")
 
# Reproject (transforms coordinates)
gdf_projected = gdf.to_crs("EPSG:3857")  # Web Mercator
gdf_projected = gdf.to_crs("EPSG:32633")  # UTM zone 33N
 
# Common CRS codes:
# EPSG:4326 - WGS84 (lat/lon)
# EPSG:3857 - Web Mercator
# EPSG:326XX - UTM zones

# Buffer (expand/shrink geometries)
buffered = gdf.geometry.buffer(100)  # 100 units buffer
 
# Centroid
centroids = gdf.geometry.centroid
 
# Simplify (reduce vertices)
simplified = gdf.geometry.simplify(tolerance=5, preserve_topology=True)
 
# Convex hull
hull = gdf.geometry.convex_hull
 
# Boundary
boundary = gdf.geometry.boundary
 
# Area and length (use projected CRS!)
gdf['area'] = gdf.geometry.area
gdf['length'] = gdf.geometry.length

# Join based on spatial relationship
joined = gpd.sjoin(gdf1, gdf2, predicate='intersects')
joined = gpd.sjoin(gdf1, gdf2, predicate='within')
joined = gpd.sjoin(gdf1, gdf2, predicate='contains')
 
# Nearest neighbor join
nearest = gpd.sjoin_nearest(gdf1, gdf2, max_distance=1000)

# Intersection
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
 
# Union
union = gpd.overlay(gdf1, gdf2, how='union')
 
# Difference
difference = gpd.overlay(gdf1, gdf2, how='difference')

# Merge geometries by attribute
dissolved = gdf.dissolve(by='region', aggfunc='sum')

# Clip data to boundary
clipped = gpd.clip(gdf, boundary_gdf)

import matplotlib.pyplot as plt
 
# Basic plot
gdf.plot()
 
# Choropleth map
gdf.plot(column='population', cmap='YlOrRd', legend=True)
 
# Multi-layer map
fig, ax = plt.subplots(figsize=(10, 10))
gdf1.plot(ax=ax, color='blue', alpha=0.5)
gdf2.plot(ax=ax, color='red', alpha=0.5)
plt.savefig('map.png', dpi=300, bbox_inches='tight')
 
# Interactive map (requires folium)
gdf.explore(column='population', legend=True)

# Join points to polygons
points_in_polygons = gpd.sjoin(points_gdf, polygons_gdf, predicate='within')
 
# Aggregate by polygon
aggregated = points_in_polygons.groupby('index_right').agg({
    'value': 'sum',
    'count': 'size'
})
 
# Merge back to polygons
result = polygons_gdf.merge(aggregated, left_index=True, right_index=True)

# Create buffers around points
gdf_projected = points_gdf.to_crs("EPSG:3857")  # Project first!
gdf_projected['buffer'] = gdf_projected.geometry.buffer(1000)  # 1km buffer
gdf_projected = gdf_projected.set_geometry('buffer')
 
# Find features within buffer
within_buffer = gpd.sjoin(other_gdf, gdf_projected, predicate='within')