ch06 corrections - choosing crs

06-reproj.qmd

@@ -18,6 +18,7 @@ import rasterio.warp
 from rasterio.plot import show
 import pyproj
 import shutil
+import math
 ```
 
 ```{python}
@@ -328,58 +329,61 @@ The [**spherely**](https://github.com/benbovy/spherely) package, which is in ear
 
 ## When to reproject? {#sec-when-to-reproject}
 
-The previous section showed how to set the CRS manually, with `lnd_layer.set_crs(4326)`. 
+The previous section showed how to set the CRS manually, with an expression such as `lnd_layer.set_crs(4326)`. 
 In real world applications, however, CRSs are usually set automatically when data is read-in. 
-In many projects the main CRS-related task is to transform objects, from one CRS into another. But when should data be transformed? And into which CRS? There are no clear-cut answers to these questions and CRS selection always involves trade-offs (Maling 1992, add reference...). 
+In many projects the main CRS-related task is to transform objects, from one CRS into another. 
+But when should data be transformed? 
+And into which CRS? 
+There are no clear-cut answers to these questions and CRS selection always involves trade-offs [@maling_coordinate_1992]. 
 However, there are some general principles provided in this section that can help you decide.
 
-First it's worth considering when to transform. In some cases transformation to a geographic CRS is essential, such as when publishing data online (for example, a Leaflet-based map using Python package `folium`). 
+First, it's worth considering when to transform. 
+In some cases transformation to a geographic CRS is essential, such as when publishing data online (for example, a Leaflet-based map using Python package [**folium**](https://python-visualization.github.io/folium/latest/)). 
 Another case is when two objects with different CRSs must be compared or combined, as shown when we try to find the distance between two objects with different CRSs:
 
 ```{python}
 lnd_layer.distance(lnd_layer_proj)
 ```
 
-We got a meaningless result and a warning.
+We got a meaningless result, and a warning.
 
-To make the `lnd_layer` and `lnd_layer_proj` objects geographically comparable one of them must be transformed into the CRS of the other. 
+To make the `lnd_layer` and `lnd_layer_proj` objects geographically comparable, one of them must be transformed into the CRS of the other. 
 But which CRS to use? 
 The answer depends on context: many projects, especially those involving web mapping, require outputs in `EPSG:4326`, in which case it is worth transforming the projected object. 
 If, however, the project requires geometric calculations, implying planar geometry, e.g., to calculating buffers (@sec-geometry-operations-on-projected-and-unprojected-data), it is necessary to transform data with a geographic CRS into an equivalent object with a projected CRS, such as the British National Grid (`EPSG:27700`). 
 That is the subject of @sec-reprojecting-raster-geometries.
 
 ## Which CRS to use? {#sec-which-crs-to-use}
 
-The question of which CRS is tricky, and there is rarely a "right" answer: "There exist no all-purpose projections, all involve distortion when far from the center of the specified frame" (Bivand, Pebesma, and Gómez-Rubio 2013, add citation...). 
+The question of which CRS is tricky, and there is rarely a "right" answer: "There exist no all-purpose projections, all involve distortion when far from the center of the specified frame" [@bivand_applied_2013]. 
 Additionally, you should not be attached just to one projection for every task. 
 It is possible to use one projection for some part of the analysis, another projection for a different part, and even some other for visualization. 
 Always try to pick the CRS that serves your goal best!
 
-When selecting **geographic** CRSs, the answer is often [WGS84](https://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS_84). 
+When selecting *geographic* CRSs, the answer is often [WGS84](https://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS_84). 
 It is used not only for web mapping, but also because GPS datasets and thousands of raster and vector datasets are provided in this CRS by default. 
-WGS84 is the most common CRS in the world, so it is worth knowing its EPSG code: 4326. 
+WGS84 is the most common CRS in the world, so it is worth knowing its EPSG code: `4326`. 
 This "magic number" can be used to convert objects with unusual projected CRSs into something that is widely understood.
 
-What about when a **projected** CRS is required? 
-In some cases, it is not something that we are free to decide: "often the choice of projection is made by a public mapping agency" (Bivand, Pebesma, and Gómez-Rubio 2013, add citation...). 
+What about when a *projected* CRS is required? 
+In some cases, it is not something that we are free to decide: "often the choice of projection is made by a public mapping agency" [@bivand_applied_2013]. 
 This means that when working with local data sources, it is likely preferable to work with the CRS in which the data was provided, to ensure compatibility, even if the official CRS is not the most accurate. 
 The example of London was easy to answer because: 
 
-* the British National Grid (with its associated EPSG code 27700) is well known, and 
+* the British National Grid (with its associated EPSG code `27700`) is well known, and 
 * the original dataset (`lnd_layer`) already had that CRS.
 
 A commonly used default is Universal Transverse Mercator ([UTM](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system)), a set of CRSs that divides the Earth into 60 longitudinal wedges and 20 latitudinal segments. 
 The transverse Mercator projection used by UTM CRSs is conformal but distorts areas and distances with increasing severity with distance from the center of the UTM zone. 
 Documentation from the GIS software Manifold therefore suggests restricting the longitudinal extent of projects using UTM zones to 6 degrees from the central meridian (source: [manifold.net](http://www.manifold.net/doc/mfd9/universal_transverse_mercator_projection.htm)). 
-Therefore, we recommend using UTM only when your focus is on preserving angles for relatively small area!
+Therefore, we recommend using UTM only when your focus is on preserving angles for a relatively small area!
 
-Almost every place on Earth has a UTM code, such as "60H" which refers to northern New Zealand. 
-UTM EPSG codes run sequentially from 32601 to 32660 for northern hemisphere locations and from 32701 to 32760 for southern hemisphere locations.
+Almost every place on Earth has a UTM code, such as `'60H'` which refers to northern New Zealand. 
+UTM EPSG codes run sequentially from `32601` to `32660` for northern hemisphere locations and from `32701` to `32760` for southern hemisphere locations.
 
-To show how the system works, let's create a function, `lonlat2UTM` to calculate the EPSG code associated with any point on the planet as follows:
+To show how the system works, let's create a function, `lonlat2UTM` to calculate the EPSG code associated with any point on the planet:
 
 ```{python}
-import math
 def lonlat2UTM(lon, lat):
  utm = (math.floor((lon + 180) / 6) % 60) + 1
  if lat > 0:
@@ -408,7 +412,7 @@ Important note: while these tools are helpful in many situations, you need to be
 In cases where an appropriate CRS is not immediately clear, the choice of CRS should depend on the properties that are most important to preserve in the subsequent maps and analysis. 
 All CRSs are either equal-area, equidistant, conformal (with shapes remaining unchanged), or some combination of compromises of those (@sec-projected-coordinate-reference-systems). 
 Custom CRSs with local parameters can be created for a region of interest and multiple CRSs can be used in projects when no single CRS suits all tasks. 
-"Geodesic calculations" can provide a fall-back if no CRSs are appropriate (see proj.org/geodesic.html). 
+"Geodesic calculations" can provide a fall-back if no CRSs are appropriate (see <https://proj.org/geodesic.html>). 
 Regardless of the projected CRS used, the results may not be accurate for geometries covering hundreds of kilometers.
 
 When deciding on a custom CRS, we recommend the following:
@@ -428,7 +432,7 @@ Note that this approach should be used with caution: no other datasets will be c
 The principles outlined in this section apply equally to vector and raster datasets. 
 Some features of CRS transformation however are unique to each geographic data model. 
 We will cover the particularities of vector data transformation in @sec-reprojecting-vector-geometries and those of raster transformation in @sec-reprojecting-raster-geometries. 
-Next, the last section, shows how to create custom map projections (@sec-custom-map-projections).
+The last section, shows how to create custom map projections (@sec-custom-map-projections).
 
 ## Reprojecting vector geometries {#sec-reprojecting-vector-geometries}