diff --git a/07-reproj.qmd b/07-reproj.qmd index c4c82c68..f94be6ac 100644 --- a/07-reproj.qmd +++ b/07-reproj.qmd @@ -674,7 +674,86 @@ axes[1].set_title('Reprojected (EPSG:32612)'); ## Custom map projections {#sec-custom-map-projections} -... +Established CRSs captured by `AUTHORITY:CODE` identifiers such as `EPSG:4326` are well suited for many applications. +However, it is desirable to use alternative projections or to create custom CRSs in some cases. +@sec-which-crs-to-use mentioned reasons for using custom CRSs, and provided several possible approaches. +Here, we show how to apply these ideas in Python. + +One is to take an existing WKT definition of a CRS, modify some of its elements, and then use the new definition for reprojecting, using the reprojection methods shown above for vector layers (@sec-reprojecting-vector-geometries) and rasters (@sec-reprojecting-raster-geometries). + +For example, let's transforms the `zion.gpkg` vector layer to a custom azimuthal equidistant (AEQD) CRS. +Using a custom AEQD CRS requires knowing the coordinates of the center point of a dataset in degrees (geographic CRS). +In our case, this information can be extracted by calculating a centroid of the `zion` area and transforming it into WGS84: + +```{python} +zion_centr = zion.centroid +zion_centr_wgs84 = zion_centr.to_crs(4326) +coords = list(zion_centr_wgs84.iloc[0].coords) +coords +``` + +Next, we can use the newly obtained lon/lat coordinates to update the WKT definition of the azimuthal equidistant (AEQD) CRS seen below. +Notice that we modified just two values below---`"Central_Meridian"` to the longitude and `"Latitude_Of_Origin"` to the latitude of our centroid: + +```{python} +my_wkt = '''PROJCS["Custom_AEQD", + GEOGCS["GCS_WGS_1984", + DATUM["WGS_1984", + SPHEROID["WGS_1984",6378137.0,298.257223563]], + PRIMEM["Greenwich",0.0], + UNIT["Degree",0.0174532925199433]], + PROJECTION["Azimuthal_Equidistant"], + PARAMETER["Central_Meridian",-113.0263], + PARAMETER["Latitude_Of_Origin",37.29818], + UNIT["Meter",1.0]]''' +``` + +This approach's last step is to transform our original object (`zion`) to our new custom CRS (`zion_aeqd`): + +```{python} +zion_aeqd = zion.to_crs(my_wkt) +``` + +Custom projections can also be made interactively, for example, using the [Projection Wizard](https://projectionwizard.org/#) web application (Šavrič, Jenny, and Jenny 2016 to add citation...). +This website allows you to select a spatial extent of your data and a distortion property, and returns a list of possible projections. +The list also contains WKT definitions of the projections that you can copy and use for reprojections. +See Open Geospatial Consortium ([2019](https://r.geocompx.org/references.html#ref-opengeospatialconsortium_wellknown_2019)) for details on creating custom CRS definitions with WKT strings. + +PROJ strings can also be used to create custom projections, accepting the limitations inherent to projections, especially of geometries covering large geographic areas, mentioned in @sec-coordinate-reference-systems. +Many projections have been developed and can be set with the +proj= element of PROJ strings, with dozens of projects described in detail on the [PROJ website](https://proj.org/operations/projections/index.html) alone. + +When mapping the world while preserving area relationships the Mollweide projection, illustrated in @fig-mollweide, is a popular and often sensible choice (Jenny et al. 2017 to add citation...). +To use this projection, we need to specify it using the proj-string element, `'+proj=moll'`, in the `.to_crs` method: + +```{python} +#| label: fig-mollweide +#| fig-cap: Mollweide projection of the world + +world.to_crs('+proj=moll').plot(color='none', edgecolor='black'); +``` + +It is often desirable to minimize distortion for all spatial properties (area, direction, distance) when mapping the world. +One of the most popular projections to achieve this is [Winkel tripel](http://www.winkel.org/other/Winkel%20Tripel%20Projections.htm), illustrated in @fig-wintri: + +```{python} +#| label: fig-wintri +#| fig-cap: Winkel tripel projection of the world + +world.to_crs('+proj=wintri').plot(color='none', edgecolor='black'); +``` + +Moreover, proj-string parameters can be modified in most CRS definitions, for example the center of the projection can be adjusted using the `+lon_0` and `+lat_0` parameters. +The below code transforms the coordinates to the Lambert azimuthal equal-area projection centered on the longitude and latitude of New York City (@fig-azimuthal-equal-area). + +```{python} +#| label: fig-azimuthal-equal-area +#| fig-cap: Lambert azimuthal equal-area projection of the world centered on New York City + +world.to_crs('+proj=laea +x_0=0 +y_0=0 +lon_0=-74 +lat_0=40') \ + .plot(color='none', edgecolor='black'); +``` + +More information on CRS modifications can be found in the [Using PROJ](https://proj.org/usage/index.html) documentation. ## Exercises