From 634af49010f1df39e81cafd14d0312dbbee08f8c Mon Sep 17 00:00:00 2001 From: Michael Dorman Date: Mon, 28 Aug 2023 23:47:06 +0200 Subject: [PATCH 1/2] multipanel and adding figure captions --- 02-spatial-data.qmd | 43 +++++++++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/02-spatial-data.qmd b/02-spatial-data.qmd index d1f51aa1..3550f814 100644 --- a/02-spatial-data.qmd +++ b/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. From 3f641738c1f78432a8e91605e8eb4521a5b892bb Mon Sep 17 00:00:00 2001 From: Michael Dorman Date: Tue, 29 Aug 2023 00:18:50 +0200 Subject: [PATCH 2/2] multipanel plot --- 04-spatial-operations.qmd | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/04-spatial-operations.qmd b/04-spatial-operations.qmd index 96023cc8..c63e0110 100644 --- a/04-spatial-operations.qmd +++ b/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