GeoPandas First Steps: Points, Coordinate Systems, Measurements and Maps

import pandas as pd
import geopandas as gpd

# Step 1: an ordinary table of cities. lat = north/south position, lon = east/west position.
cities = pd.DataFrame({
    "city": ["Bengaluru", "Mysuru", "Mangaluru", "Hubballi", "Chennai", "Hyderabad", "Mumbai", "Delhi"],
    "state": ["Karnataka", "Karnataka", "Karnataka", "Karnataka", "Tamil Nadu", "Telangana", "Maharashtra", "Delhi"],
    "population": [8443675, 920550, 484785, 943857, 4646732, 6809970, 12442373, 11034555],
    "lat": [12.9716, 12.2958, 12.9141, 15.3647, 13.0827, 17.3850, 19.0760, 28.6139],
    "lon": [77.5946, 76.6394, 74.8560, 75.1240, 80.2707, 78.4867, 72.8777, 77.2090],
})

# Step 2: turn it into a map layer.
# points_from_xy makes a point from each pair. The order is always x first (lon), then y (lat).
# crs="EPSG:4326" says the numbers are longitude and latitude (GPS style coordinates).
gdf = gpd.GeoDataFrame(
    cities,
    geometry=gpd.points_from_xy(cities["lon"], cities["lat"]),
    crs="EPSG:4326",
)
print(gdf[["city", "geometry"]].head(3))      # a new "geometry" column holds the points
print(type(gdf).__name__)

        city                 geometry
0  Bengaluru  POINT (77.5946 12.9716)
1     Mysuru  POINT (76.6394 12.2958)
2  Mangaluru   POINT (74.856 12.9141)
GeoDataFrame

print(gdf.crs)                                    # the coordinate system of the layer
print(gdf.geom_type.unique().tolist())            # the kinds of shape in it (here only Point)
print(gdf.total_bounds.round(2))                  # the box around everything: min x, min y, max x, max y
print(gdf.geometry.iloc[0])                       # the first shape, in text form
print(gdf.geometry.iloc[0].x, gdf.geometry.iloc[0].y)     # its x and y on their own

EPSG:4326
['Point']
[72.88 12.3  80.27 28.61]
POINT (77.5946 12.9716)
77.5946 12.9716

# Step 1: keep only Karnataka's cities. .copy() makes a separate table we can change safely.
karnataka = gdf[gdf["state"] == "Karnataka"].copy()

# Step 2: to_crs() converts the coordinates to another system.
# EPSG:32643 is UTM zone 43N, which measures in metres (good for distance and area).
utm = karnataka.to_crs("EPSG:32643")

# Same place, two sets of numbers: degrees before, metres after
print(karnataka.crs.to_epsg(), utm.crs.to_epsg())
print(karnataka.geometry.iloc[0])
print(utm.geometry.iloc[0])

4326 32643
POINT (77.5946 12.9716)
POINT (781481.0202947541 1435426.4069104537)

# Step 1: pick out Bengaluru's point
blr = utm[utm["city"] == "Bengaluru"].geometry.iloc[0]

# Step 2: distance from Bengaluru to every city. In UTM the answer is in metres, so divide by 1000 for km.
utm["km_from_bengaluru"] = (utm.geometry.distance(blr) / 1000).round(1)
print(utm[["city", "km_from_bengaluru"]])

        city  km_from_bengaluru
0  Bengaluru                0.0
1     Mysuru              127.9
2  Mangaluru              297.2
3   Hubballi              375.8

from pyproj import Geod                           # a tool that measures on the curved Earth

geod = Geod(ellps="WGS84")                        # the Earth shape used by GPS

# Two places given as longitude and latitude: Bengaluru and Mysuru
lon1, lat1 = 77.5946, 12.9716
lon2, lat2 = 76.6394, 12.2958

# inv() returns three values: direction there, direction back, distance in metres.
# We only want the distance, so the first two are thrown away with _
_, _, metres = geod.inv(lon1, lat1, lon2, lat2)
print(round(metres / 1000, 1))                    # compare with the UTM answer above

127.9

from shapely.geometry import box

# A box around Bengaluru, written in degrees: box(west, south, east, north)
area = gpd.GeoDataFrame(
    {"name": ["test box"]},
    geometry=[box(77.4, 12.8, 77.8, 13.2)],
    crs="EPSG:4326",
)

# In degrees the "area" is a meaningless number (degrees squared)
print("in degrees squared:", round(area.geometry.area.iloc[0], 4))

# Convert to metres first, then the area and length are real
area_utm = area.to_crs("EPSG:32643")
print("in km2:", round(area_utm.geometry.area.iloc[0] / 1_000_000, 1))
print("perimeter km:", round(area_utm.geometry.length.iloc[0] / 1000, 1))

in degrees squared: 0.16
in km2: 1922.5
perimeter km: 175.4

import matplotlib.pyplot as plt                   # the drawing library

# Step 1: create an empty picture (fig) with one drawing area (ax)
fig, ax = plt.subplots(figsize=(9, 8))

# Step 2: draw the points. Colour and size both follow the population.
gdf.plot(
    ax=ax,
    column="population",                          # which column decides the colour
    cmap="YlOrRd",                                # colour ramp: yellow to red
    markersize=gdf["population"] / 40_000,        # bigger city, bigger dot
    edgecolor="black",
    legend=True,
    legend_kwds={"label": "Population", "shrink": 0.6},
)
ax.margins(0.12)                                  # add some empty space around the points

# Step 3: write each city name next to its dot
for _, row in gdf.iterrows():                     # go through the table one row at a time
    ax.annotate(row["city"], (row["lon"], row["lat"]), xytext=(6, 5), textcoords="offset points", fontsize=9)

# Step 4: title, axis names, and save to a file
ax.set_title("Cities by population")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
fig.savefig("cities_map.png", dpi=150, bbox_inches="tight")
print("saved")

saved

# Save the layer in two common GIS formats. GeoPackage (.gpkg) is the better default.
gdf.to_file("cities.gpkg", layer="cities", driver="GPKG")
gdf.to_file("cities.geojson", driver="GeoJSON")

# Read it back to prove the save worked: same number of rows, same coordinate system
back = gpd.read_file("cities.gpkg", layer="cities")
print(len(back), back.crs.to_epsg())

8 4326