spatialdataPy.qmd

---
format:
  html:
    code-fold: false
jupyter: python3
---

# Lab in Python {#sec-spatial-data-Python .unnumbered}

In this lab, we will learn how to load, manipulate and visualize spatial data. In some senses, spatial data are usually included simply as "one more column" in a table. However, *spatial is special* sometimes and there are few aspects in which geographic data differ from standard numerical tables. In this session, we will extend the skills developed in the previous one about non-spatial data, and combine them. In the process, we will discover that, although with some particularities, dealing with spatial data in `Python` largely resembles dealing with non-spatial data.

## Importing modules

```{python}
#| warning: false

# Import the pandas library, which is useful for data manipulation and analysis.
import pandas as pd
# Import the geopandas library, which extends pandas to support spatial data operations.
import geopandas as gpd
# Import the Point class from the shapely.geometry module, used for handling geometric points.
from shapely.geometry import Point
# Import the osmnx library, which simplifies the process of downloading and analyzing street networks and other geospatial data from OpenStreetMap.
import osmnx as ox
# Import the contextily library, which is used for adding basemaps (like raster tiles) to geospatial plots.
import contextily as cx
# Import the pyplot module from matplotlib, a library for creating static, animated, and interactive visualizations in Python.
import matplotlib.pyplot as plt
# Import the CRS class from pyproj, which provides tools for defining and transforming coordinate reference systems.
from pyproj import CRS
# Operating systems
import os
# Interactive maps
import folium
# Import the display function
from IPython.display import display

```

## Datasets

Today we are going to go to London. We will be playing around with different datasets loading them both locally and dynamically from the web. You can download data manually, keep a copy on your computer, and load them from there.

### Creating geographic data

First we will use the following commands create geographic datasets *from scratch* representing coordinates of some famous locations in London. Most projects start with pre-generated data, but it's useful to create datasets to understand data structures.

```{python}
# Create the DataFrame
data = {
    'name': ["The British Museum", "Big Ben", "King's Cross", "The Natural History Museum"],
    'lon': [-0.1459604, -0.1272057, -0.1319481, -0.173734],
    'lat': [51.5045975, 51.5007325, 51.5301701, 51.4938451]
}
poi_df = pd.DataFrame(data)

# Convert DataFrame to GeoDataFrame
geometry = [Point(xy) for xy in zip(poi_df['lon'], poi_df['lat'])]
poi_gdf = gpd.GeoDataFrame(poi_df, geometry=geometry)

# Set the coordinate reference system (CRS)
poi_gdf.set_crs(epsg=4326, inplace=True)

print(poi_gdf)
```

### Types of Data

Now let's look at the different types of geographical data starting with polygons. We will use a dataset that contains the boundaries of the districts of London. We can read it into an object named districts.

::: {.panel-tabset group="data"}
## Polygons

We first import the district shapefile use `gpd.read_file`, we then plot it to make sure we are seeing it 'correctly'.

```{python}
# Read the shapefile for districts
districts = gpd.read_file("data/London/Polygons/districts.shp")

# Create a simple plot
districts.plot()

# Display the plot
plt.show()
```

## Lines

We them import a file of roads in London and plot it.

```{python}
# Read the shapefile for A roads
a_roads = gpd.read_file("data/London/Lines/a_roads.shp")

# If you needed to import a `geojson` file, this would be the function:
# a_roads = gpd.read_file("data/London/Lines/a_roads.geojson")

# Create a simple plot of the roads
a_roads.plot()

# Display the plot
plt.show()
```

## Points

We can also import point files. So far, we have imported `shapefiles` and `geojsons`, but we can also obtain data from urls like in the [Open Science DIY](https://pietrostefani.github.io/gds/openscienceDIY.html) session or from other sources like **OpenStreetMap**. Both `R` and `Python` have libraries that allow us to query OpenStreetMap.

Note that we use the method `features_from_place`, which queries for points in a particular place (London in this case) and creates a GeoDataFrame of OSM features.

```{python}
# Create an OSM query for "Greater London, U.K."
query = "London, United Kingdom"
restaurants = ox.features_from_place(query, tags={"amenity": ["restaurant", "bar", "pub"]})
# Create a simple plot of the roads
restaurants.plot()
# Display the plot
plt.show()

```

And to inspect the data queried:

```{python}
restaurants.info()
```

You do not need to know at this point what happens behind the scenes when we run these lines but, if you are curious, we are making a query to OpenStreetMap (almost as if you typed "restaurant in London, UK" within Google Maps) and getting the response as a table of data, instead of as a website with an interactive map. Pretty cool, huh?

*Note*: the code cells above requires internet connectivity. For more about querying from osm see [here](https://osmnx.readthedocs.io/en/stable/user-reference.html#osmnx.features.features_from_point).

**Important**: Be careful, if you query too much data, your environment is likely to get stuck.
:::

## Inspecting Spatial Data

### Inspecting

Just like a `dataframe` (see the OpenScience Lab), we can inspect the data (attributes table) within a spatial object. The most direct way to get from a file to a quick visualization of the data is by loading it and calling the `plot` command. Let's start by inspecting the data like we did for non spatial `dataframes`.

We can see our data is very similar to a traditional, non-spatial `dataFrame`, but with an additional column called geometry.

```{python}
# Read the first 5 rows of the data
print(districts.head()) 
```

We can inspect the object in different ways :

```{python}
# Read the first row
print(districts.iloc[0])

# Read the first column
print(districts.iloc[:, 0])

# Read the first row, first column
print(districts.iloc[0, 0])

# Read the column "DIST_NAME"
print(districts['DIST_NAME'])
```

We can read or create subsets:

```{python}
# dataframe can be subsetted using conditional statement
# read the rows which have "City of London" as value for DIST_NAME
# Filter rows where 'DIST_NAME' is 'City of London'
filtered_districts = districts[districts['DIST_NAME'] == 'City of London']

print(filtered_districts)
```

### Quick visualisation

Let's start by plotting London in a colour and adding Hackney (a district) in a different colour.

```{python}
# Plot London in grey
fig, ax = plt.subplots()
districts.plot(ax=ax, color='lightgrey')

# Add city of London (Hackney) in turquoise to the map
hackney = districts[districts['DIST_NAME'] == 'Hackney']
hackney.plot(ax=ax, color='turquoise')

plt.show()
```

Some guidance on colours in `Python` can be found [here](https://matplotlib.org/stable/gallery/color/named_colors.html).

## Styling plots

It is possible to tweak many aspects of a plot to customize if to particular needs. In this section, we will explore some of the basic elements that will allow us to obtain more compelling maps.

**Note**: some of these variations are very straightforward while others are more intricate and require tinkering with the internal parts of a plot. They are not necessarily organized by increasing level of complexity.

### Plotting different layers

We first start by plotting one layer over another

```{python}
# Plotting the geometries
fig, ax = plt.subplots()

# Plot districts with no fill (transparent fill)
districts.plot(ax=ax, edgecolor='black', facecolor='none')  # No fill, only border
# Plot roads with transparency
a_roads.plot(ax=ax, color='brown')  # Roads in brown

plt.show()
```

### Changing transparency

The intensity of color of a polygon can be easily changed through the alpha attribute in plot. This is specified as a value betwee zero and one, where the former is entirely transparent while the latter is the fully opaque (maximum intensity):

```{python}
# Plotting the geometries
fig, ax = plt.subplots()

# Plot districts with no fill (transparent fill)
districts.plot(ax=ax, edgecolor='black', facecolor='none')  # No fill, only border
# Plot roads with transparency
a_roads.plot(ax=ax, color='brown', alpha=0.5)  # Roads in brown

plt.show()
```

### Removing axes

Although in some cases, the axes can be useful to obtain context, most of the times maps look and feel better without them. Removing the axes involves wrapping the plot into a figure, which takes a few more lines of aparently useless code but that, in time, it will allow you to tweak the map further and to create much more flexible designs.

```{python}
# Plotting the geometries
fig, ax = plt.subplots()

# Plot districts with no fill (transparent fill)
districts.plot(ax=ax, edgecolor='black', facecolor='none')  # No fill, only border
# Plot roads with transparency
a_roads.plot(ax=ax, color='brown', alpha=0.5)  # Roads with 50% transparency

# Remove the axis
ax.axis('off')

plt.show()
```

Let us stop for a second a study each of the previous lines:

We have first created a figure named `fig` with one axis named `ax` by using the command `plt.subplots` (part of the library matplotlib, which we have imported at the top of the notebook). Note how the method is returning two elements and we can assign each of them to objects with different name (`fig` and `ax`) by simply listing them at the front of the line, separated by commas.

Second, we plot the geographies as before, but this time we tell the function that we want it to draw the polygons on the axis we are passing, `ax`. This method returns the axis with the geographies in them, so we make sure to store it on an object with the same name, `ax`.

On the third line, we effectively remove the box with coordinates.

Finally, we draw the entire plot by calling plt.`show()`.

### Adding a title

Adding a title is an extra line, if we are creating the plot within a figure, as we just did. To include text on top of the figure:

```{python}
# Plotting the geometries
fig, ax = plt.subplots()

# Plot districts with no fill (transparent fill)
districts.plot(ax=ax, edgecolor='black', facecolor='none')  # No fill, only border
# Plot roads with transparency
a_roads.plot(ax=ax, color='brown', alpha=0.5)  # Roads with 50% transparency
# Remove the axis
ax.axis('off')
# Add figure title
fig.suptitle("Main roads in London")
# Display
plt.show()
```

### Changing what border lines look like

Border lines sometimes can distort or impede proper interpretation of a map. In those cases, it is useful to know how they can be modified. Let us first see the code to make the lines thicker and black, and then we will work our way through the different steps:

```{python}
# Convert CRS from WGS 84 (EPSG:4326) to British National Grid (EPSG:27700)
poi_gdf_bng = poi_gdf.to_crs(epsg=27700)

# Plotting the geometries
fig, ax = plt.subplots()
# Plot districts with no fill, black borders
districts.plot(ax=ax, edgecolor='black', facecolor='none')
# Plot roads with brown color and 50% transparency
a_roads.plot(ax=ax, color='brown', alpha=0.5)
# Plot restaurants with blue color and adjusted size
poi_gdf_bng.plot(ax=ax, edgecolor='blue', facecolor='blue', markersize=100)# Adjust size accordingly
# Remove the axis for a clean look
ax.axis('off')
# Add figure title
fig.suptitle("Main roads in London")
# Display the plot
plt.show()

```

### Labelling

Labeling maps is of paramount importance as it is often key when presenting data analysis and visualization. Properly labeled maps enables readers to effectively analyze and interpret spatial data.

-   `iterrows()` is a pandas function that iterates over DataFrame rows as (index, Series) pairs. Here, `idx` is the index of the row, and row is a `pandas` series containing the data for that particular district.
-   `centroid = row['geometry'].centroid` gets the centroid of each district's geometry. `row['geometry']` refers to the geometry of the district, which could be a polygon or a multipolygon. `.centroid` computes the geometric center (centroid) of this polygon.
-   `ax.text()` is a method from `Matplotlib`, used to place text at specific coordinates on the plot. `centroid.x` and `centroid.y` provide the x and y coordinates of the centroid, which determine where the text will be placed. `row['DIST_NAME']` is the name of the district that will be displayed as the label. `fontsize=6` sets the size of the text to 8 points. `ha='center'` ensures that the text is horizontally aligned to the center of the specified coordinates.

```{python}
# Plot the districts with a gray fill
fig, ax = plt.subplots()
districts.plot(ax=ax, edgecolor="black", facecolor='none')

# Add text labels at the centroids of the districts
for idx, row in districts.iterrows():
    centroid = row['geometry'].centroid
    ax.text(centroid.x, centroid.y, row['DIST_NAME'], fontsize=6, ha='center')

# Remove axis
ax.set_axis_off()

plt.show()

```

## Coordinate reference Systems

### CRSs in Python

Coordindate reference systems (CRS) are the way geographers and cartographers represent a three-dimentional objects, such as the round earth, on a two-dimensional plane, such as a piece of paper or a computer screen. If the source data contain information on the CRS of the data, we can modify this.

First we need to retrieve the CRS from the vector data.

```{python}
# Retrieve the CRS from the GeoDataFrame
crs = districts.crs
# Print the CRS information
print(crs)
```

We can also retrieve some additional information about the used CRS. For example, try to run:

```{python}
# Check if the CRS is geographic or not
is_geographic = CRS(crs).is_geographic
# Find out the CRS units
units_gdal = CRS(crs).axis_info[0].unit_name if CRS(crs).axis_info else None
# Extract the SRID
srid = CRS(crs).to_epsg()
# Extract the proj4string representation
proj4string = CRS(crs).to_proj4()

# Print results
print(f"Is Geographic: {is_geographic}")
print(f"Units (GDAL): {units_gdal}")
print(f"SRID: {srid}")
print(f"Proj4 String: {proj4string}")
```

As we can see, there is information stored about the reference system: it is using the standard British projection (British National Grid EPSG:27700), which is expressed in meters. There are also other less decipherable parameters but we do not need to worry about them right now.

If we want to modify this and "reproject" the polygons into a different CRS, the quickest way is to find the EPSG code online (epsg.io is a good one, although there are others too). For example, if we wanted to transform the dataset into lat/lon coordinates, we would use its EPSG code, 4326 (CRS's name "WGS84"):

In cases when a coordinate reference system (CRS) is missing or the wrong CRS is set, the `.to_crs` function can be used:

```{python}
# Transform the CRS to EPSG:4326
districts_4326 = districts.to_crs(epsg=4326)

# Optionally, print the new CRS to verify
print(districts_4326.crs)
```

### From coordinates to spatial objects

CRSs are also very useful if we obtain data that is in a csv, has coordinates but needs to be transformed to a `GeoDataFrame`. For example we have some London housing transactions we want to import and use.

We want to transform the .csv in a `GeoDataFrame` using the coordinates stored in columns 17 and 18, and then we set the `GeoDataFrame` CRS to the British National Grid (EPSG:27700).

```{python}
# Import housesales data from CSV
housesales = pd.read_csv("data/London/Tables/housesales.csv")

# Filter housesales to include only those with price less than 500000
housesales_f = housesales[housesales['price'] < 500000]

# Assume columns 17 and 18 are 'longitude' and 'latitude' respectively
housesales_gdf = gpd.GeoDataFrame(
    housesales_f, geometry=gpd.points_from_xy(housesales_f.greastings, housesales_f.grnorthing), crs="EPSG:27700"
)

print(housesales_gdf.head())
```

### Zooming in or out

It's important to know what CRS your data is in if you want to create zoomed versions of your maps. [BBox finder](http://bboxfinder.com/#0.000000,0.000000,0.000000,0.000000) is a useful tool to identify coordinates in `EPSG:4326`.

Here for example we are zooming in to some of the point we created at the beginning of the lab.

```{python}
# Create a plot
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the districts
districts_4326.plot(ax=ax, color='none', edgecolor='black')
# Plot the points of interest
poi_gdf.plot(ax=ax, color='blue', markersize=50)
# Set the coordinate limits
ax.set_xlim(-0.180723, -0.014212)
ax.set_ylim(51.476668, 51.532337)
# Remove axis labels and ticks
ax.set_axis_off()
# Show the plot
plt.show()

```

## Manipulating Spatial Tables

Once we have an understanding of how to visually display spatial information contained, let us see how it can be combined with the operations related to manipulating non-spatial tabular data. Essentially, the key is to realize that a `GeoDataFrame` contain most of its spatial information in a single column named geometry, but the rest of it looks and behaves exactly like a non-spatial `DataFrame` (in fact, it is). This concedes them all the flexibility and convenience that we saw in manipulating, slicing, and transforming tabular data, with the bonus that spatial data is carried away in all those steps. In addition, `GeoDataFrame` also incorporate a set of explicitly spatial operations to combine and transform data. In this section, we will consider both.

`GeoDataFrames` come with a whole range of traditional GIS operations built-in. Here we will run through a small subset of them that contains some of the most commonly used ones.

::: {.panel-tabset group="data"}
## Area

One of the spatial aspects we often need from polygons is their area. "How big is it?" is a question that always haunts us when we think of countries, regions, or cities. To obtain area measurements, first make sure the `GeoDataFrame` you are working with is projected. If that is the case, you can calculate areas as follows:

We had already checked that district was projected to the British National Grid

```{python}
districts_areas = districts.area
districts_areas.head()

areas_in_sqkm = districts_areas / 1000000 #convert into squared kilometres
areas_in_sqkm.head()
```

## Length

Similarly, an equally common question with lines is their length. Also similarly, their computation is relatively straightforward, provided that our data are projected.

```{python}
street_length = a_roads.length
street_length.head()
```

If you check the `dataframe` you will see the lengths.

## Centroids

Sometimes it is useful to summarize a polygon into a single point and, for that, a good candidate is its centroid (almost like a spatial analogue of the average).

```{python}
# Create a dataframe with centroids
cents = districts.centroid
cents.head()
```

```{python}
# Create a plot
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the districts
districts.plot(ax=ax, color='lightgrey', edgecolor='black')
# Plot the centroids
cents.plot(ax=ax, color='red', markersize=10, edgecolor='none')
# Set minimal theme by removing axes and grid
ax.set_axis_off()
# Show the plot
plt.show()
```

## Buffers and selecting by location

To create a buffer using `geopandas`, simply call the buffer method, passing in the radious. For example, to draw a 1000m. buffer around every centroid of every district:

```{python}
# buffer
buf = cents.buffer(1000)
buf.head()
```

```{python}
# Create a plot
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the districts
districts.plot(ax=ax, color='none', edgecolor='black')
# Plot the centroids
cents.plot(ax=ax, color='black', markersize=20, edgecolor='none')
# Plot the centroid buffers
buf.plot(ax=ax, color='none', alpha=0.5, edgecolor='red')
# Set minimal theme by removing axes and grid
ax.set_axis_off()
# Show the plot
plt.show()
```
:::

## Joins

## Join districts with educational level data

```{python}
# Import qualifications data from CSV
qualifications2001_df = pd.read_csv("data/London/Tables/qualifications2001_2.csv")

# Take a quick look at the table by reading the first 5 lines
qualifications2001_df.head()
```

Check the data

```{python}
# Join check that the code you are joining the data on is the same.
# First check the GeoDataFrame
districts.head()

# Then the DataFrame
qualifications2001_df.head()

## rename(columns={'Zone-Code': 'Dist_code'}) specifies that the column name Zone-Code should be renamed to Dist_code
qualifications2001_df.rename(columns={'Zone_Code': 'DIST_CODE'}, inplace=True)
```

Merge the qualification dataset to the districts `GeoDataFrame`:

-   Merge the `qualifications2001_df` DataFrame with the `districts` `GeoDataFrame`.
-   The merge is done on the `DIST_CODE` column, which must be present in both DataFrames.
-   By default, this performs an inner join, but you can specify different types of joins (e.g., left, right) with the 'how' parameter.
-   `gpd.GeoDataFrame(...)` converts the resulting merged DataFrame into a GeoDataFrame using the 'geopandas' library.
-   `district_qual` is the new `GeoDataFrame` that contains the combined data from 'qualifications2001_df' and 'districts'.
-   `districts.plot()` plots the `GeoDataFrame` and `column='level 4'` specifies the column to plot.
-   `legend=True` adds a legend to the plot and `cmap='OrRd'` applies a color map.

```{python}
district_qual = districts.merge(qualifications2001_df, on='DIST_CODE')
# Prove it worked
district_qual.plot(column="Level4", cmap="OrRd")
# Show the plot
plt.show()

```

### Calculation

Now, let's create the share of people with level 4 qualification, i.e. create the new variable `Level4p` equal to the number of people with level4 qualification divided by total population:

```{python}
# Assuming 'Level4' and 'Pop' columns exist in the merged GeoDataFrame district_qual
district_qual['Level4p'] = district_qual['Level4'] / district_qual['Population1674']

```

## Saving maps to figures

Create a file to put your maps:

```{python}
# Create directory "maps" if it doesn't exist
os.makedirs("maps2", exist_ok=True)
```

Let's create a simple map with the variable we just created and save to external file:

```{python}
# Create a figure and axes for the plot
fig, ax = plt.subplots()
# Plot the districts' geometry with no fill color and black edges
districts.plot(ax=ax, color='none', edgecolor='black')
# Plot the district qualifications with a color map
district_qual.plot(ax=ax, column="Level4", cmap="OrRd", legend=True)
# Save the plot to a PDF file
plt.savefig("maps/london_test2.pdf")
# Save the plot as a JPG file
plt.savefig("maps/london_test2.jpg", format='jpg', dpi=300)
# Close the plot
plt.close()
```

## Adding baselayers

Many single datasets lack context when displayed on their own. A common approach to alleviate this is to use web tiles, which are a way of quickly obtaining geographical context to present spatial data. In Python, we can use `contextily` to pull down tiles and display them along with our own geographic data.

We can begin by creating a map in the same way we would do normally, and then use the `add_basemap` command to add a basemap:

```{python}
ax = restaurants.plot(color="black")
cx.add_basemap(ax, crs=restaurants.crs, source=cx.providers.CartoDB.Voyager);

# Show the plot
plt.show()
```

Note that we need to be explicit when adding the basemap to state the coordinate reference system (crs) our data is expressed in, `contextily` will not be able to pick it up otherwise. Conversely, we could change our data’s CRS into `Pseudo-Mercator`, the native reference system for most web tiles:

```{python}
restaurants_wm = restaurants.to_crs(epsg=3857)
ax = restaurants_wm.plot(color="black")
cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik);

# Show the plot
plt.show()
```

Note how the coordinates are different but, if we set it right, either approach aligns tiles and data nicely.

Now, contextily offers a lot of options in terms of the sources and providers you can use to create your basemaps. For example, we can use satellite imagery instead. A lot more about basemap options with `contextly` [here](https://contextily.readthedocs.io/en/latest/providers_deepdive.html).

```{python}
f, ax = plt.subplots(1, figsize=(9, 9))
districts.plot(alpha=0.5, ax=ax)
cx.add_basemap(
    ax, 
    crs=districts.crs,
    source=cx.providers.Esri.WorldImagery
)
ax.set_axis_off()
# Show the plot
plt.show()
```

## Interactive maps

Everything we have seen so far relates to static maps. These are useful for publication, to include in reports or to print. However, modern web technologies afford much more flexibility to explore spatial data interactively.

We wil use the `folium` library, which is a Python wrapper for `Leaflet`. Here’s how you can do it:

```{python}
# Define the locations and popups
locations = {
    "The British Museum": (51.5045975, -0.1459604),
    "Big Ben": (51.5007325, -0.1272057),
    "King's Cross": (51.5301701, -0.1319481),
    "The Natural History Museum": (51.4938451, -0.173734)
}

# Create a base map
m = folium.Map(location=[51.505, -0.127], zoom_start=14, tiles='CartoDB positron')

# Add markers
for popup, (lat, lng) in locations.items():
    folium.Marker(location=[lat, lng], popup=popup).add_to(m)

# Display the map
display(m)
```