Merge branch 'main' of https://github.com/geocompx/geocompy

02-spatial-data.qmd

@@ -159,7 +159,7 @@ else:
 gdf = gpd.read_file('data/world.gpkg')
 ```
 
-As  result is an object of type (class) `GeoDataFrame` with 177 rows (features) and 11 columns, as shown in the output of the following code:
+As result is an object of type (class) `GeoDataFrame` with 177 rows (features) and 11 columns, as shown in the output of the following code:
 
 ```{python}
 #| label: typegdf
@@ -183,9 +183,12 @@ The following expression creates a subset based on a condition, such as equality
 gdf[gdf['name_long'] == 'Egypt']
 ```
 
-Finally, to get a sense of the spatial component of the vector layer, it can be plotted using the `.plot` method, as follows:
+Finally, to get a sense of the spatial component of the vector layer, it can be plotted using the `.plot` method (@fig-gdf-plot):
 
 ```{python}
+#| label: fig-gdf-plot
+#| fig-cap: Basic plot of a `GeoDataFrame`
+
 gdf.plot();
 ```
 
@@ -630,9 +633,12 @@ src = rasterio.open('data/srtm.tif')
 src
 ```
 
-To get a first impression of the raster values, we can plot the raster using the `show` function:
+To get a first impression of the raster values, we can plot the raster using the `show` function (@fig-rasterio-plot):
 
 ```{python}
+#| label: fig-rasterio-plot
+#| fig-cap: Basic plot of a raster, the data are coming from a `rasterio` file connection
+
 rasterio.plot.show(src);
 ```
 
@@ -740,15 +746,21 @@ new_transform
 
 Note that, confusingly, $delta_{y}$ is defined in `rasterio.transform.from_origin` using a positive value (`0.5`), even though it is eventuially negative (`-0.5`)! 
 
-The raster can now be plotted in its coordinate system, passing the array along with the transformation matrix to `show`:
+The raster can now be plotted in its coordinate system, passing the array `elev` along with the transformation matrix `new_transform` to `show` (@fig-rasterio-plot-elev):
 
 ```{python}
+#| label: fig-rasterio-plot-elev
+#| fig-cap: Plot of the `elev` raster, which was created from scratch
+
 rasterio.plot.show(elev, transform=new_transform);
 ```
 
-The `grain` raster can be plotted the same way, as we are going to use the same transformation matrix for it as well:
+The `grain` raster can be plotted the same way, as we are going to use the same transformation matrix for it as well (@fig-rasterio-plot-grain):
 
 ```{python}
+#| label: fig-rasterio-plot-grain
+#| fig-cap: Plot of the `grain` raster, which was created from scratch
+
 rasterio.plot.show(grain, transform=new_transform);
 ```
 
@@ -890,17 +902,16 @@ And here is an illustration of the layer in the original projected CRS and in a
 
 ```{python}
 #| label: fig-zion-crs
-#| fig-cap: Examples of geographic (WGS 84; left) and projected (NAD83 / UTM zone 12N; right) coordinate systems for a vector data type. 
-
-fig, axes = plt.subplots(ncols=2, figsize=(9,4))
-zion.to_crs(4326).plot(ax=axes[0], edgecolor='black', color='lightgrey')
-zion.plot(ax=axes[1], edgecolor='black', color='lightgrey')
-axes[0].set_axisbelow(True)  ## Plot grid below other elements
-axes[1].set_axisbelow(True)
-axes[0].grid()  ## Add grid
-axes[1].grid()
-axes[0].set_title('WGS 84')
-axes[1].set_title('UTM zone 12N');
+#| fig-cap: Examples of Coordinate Refrence Systems (CRS) for a vector layer 
+#| fig-subcap: 
+#| - geographic (WGS84)
+#| - projected (NAD83 / UTM zone 12N)
+#| layout-ncol: 2
+
+# WGS84
+zion.to_crs(4326).plot(edgecolor='black', color='lightgrey').grid()
+# NAD83 / UTM zone 12N
+zion.plot(edgecolor='black', color='lightgrey').grid();
 ```
 
 We are going to elaborate on reprojection from one CRS to another (`.to_crs` in the above code section) in @sec-reproj-geo-data.

04-spatial-operations.qmd

@@ -331,7 +331,7 @@ fig.delaxes(axes[1][1]);
 
 Sometimes two geographic datasets do not touch but still have a strong geographic relationship. 
 The datasets `cycle_hire` and `cycle_hire_osm`, provide a good example.
-Plotting them shows that they are often closely related but they do not touch, as shown in @fig-cycle-hire, which is created with the code below: 
+Plotting them shows that they are often closely related but they do not touch, as shown in @fig-cycle-hire:
 
 ```{python}
 #| label: fig-cycle-hire
@@ -354,24 +354,23 @@ m.to_numpy().any()
 Imagine that we need to join the capacity variable in `cycle_hire_osm` onto the official 'target' data contained in `cycle_hire`.
 This is when a non-overlapping join is needed. 
 Spatial join (`gpd.sjoin`) along with buffered geometries can be used to do that.
-This is demonstrated below, using a threshold distance of 20 m. 
+This is demonstrated below, using a threshold distance of 20 $m$. 
 Note that we transform the data to a projected CRS (`27700`) to use real buffer distances, in meters. 
 
 ```{python}
 crs = 27700
-cycle_hire2 = cycle_hire.copy().to_crs(crs)
-cycle_hire2['geometry'] = cycle_hire2.buffer(20)
+cycle_hire_buffers = cycle_hire.copy().to_crs(crs)
+cycle_hire_buffers['geometry'] = cycle_hire_buffers.buffer(20)
 z = gpd.sjoin(
-  cycle_hire2, 
+  cycle_hire_buffers, 
   cycle_hire_osm.to_crs(crs)
 )
-z[['name_left','name_right']]
+z
 ```
 
-The result `z` shows that there are 438 points in the target object `cycle_hire` within the threshold distance of `cycle_hire_osm`.
 Note that the number of rows in the joined result is greater than the target.
-This is because some cycle hire stations in `cycle_hire` have multiple matches in `cycle_hire_osm`.
-To aggregate the values for the overlapping points and return the mean, we can use the aggregation methods learned in @attr, resulting in an object with the same number of rows as the target.
+This is because some cycle hire stations in `cycle_hire_buffers` have multiple matches in `cycle_hire_osm`.
+To aggregate the values for the overlapping points and return the mean, we can use the aggregation methods learned in @sec-vector-attribute-aggregation, resulting in an object with the same number of rows as the target.
 We also go back from buffers to points using `.centroid`:
 
 ```{python}
@@ -386,13 +385,18 @@ The capacity of nearby stations can be verified by comparing a plot of the capac
 
 ```{python}
 #| label: fig-cycle-hire-z
-#| fig-cap: Non-overlapping join input (left) and result (right)
-
-fig, axes = plt.subplots(ncols=2, figsize=(8,3))
-cycle_hire_osm.plot(column='capacity', legend=True, ax=axes[0])
-z.plot(column='capacity', legend=True, ax=axes[1])
-axes[0].set_title('Input')
-axes[1].set_title('Join result');
+#| fig-cap: Non-overlapping join
+#| fig-subcap: 
+#| - Input
+#| - Join result
+#| layout-ncol: 2
+
+# Input
+fig, ax = plt.subplots(1, 1, figsize=(6, 3))
+cycle_hire_osm.plot(column='capacity', legend=True, ax=ax);
+# Join result
+fig, ax = plt.subplots(1, 1, figsize=(6, 3))
+z.plot(column='capacity', legend=True, ax=ax);
 ```
 
 ### Spatial aggregation