[Feat] Make EPSG:4087 the default CRS

[ldesousa]

Is your feature request related to a problem? Please describe. GRASS GIS currently uses the WGS:84 datum ensemble as default location, for instance in the new GUI. However, when displaying layers in this system the Equirectangular projection is applied implicitly. Moreover, the GUI reports coordinates in the wrong order for the WGS:84 (latitude, longitude).

Describe the solution you'd like GRASS GIS is a spatial analysis system operating on the Cartesian plane, most modules do not apply to the geodetic (or geographic) domain. Therefore the default location should be an explicitly projected CRS. The current behaviour, WGS:84 plus Equirectangular projection, matches the formal CRS definition in the EPSG database with code 4087.

Describe alternatives you've considered Also worth considering not using the WGS:84 ensemble, pointing instead to a concrete datum. However there is no obvious choice, and it is important to not force datum transformations importing GPS/GNSS data to the default location. Would be great if one day the EPSG could issue a code identifying the latest datum in the WGS:84 series.

To address the axes order issue exclusively, the CRS:84 system specified by the OGC is a viable alternative.

Additional context

The EPSG:4087 page at epsg.io: https://epsg.io/4087

The Wikipaedia entry on the Equirectangular projection: https://en.wikipedia.org/wiki/Equirectangular_projection

The PROJ manual on the Equirectangular projection: https://proj.org/en/9.2/operations/projections/eqc.html

[jjimenezshaw]

Hi. I was talking with @ldesousa in FOSS4G, in particular about this issue. I'm not aware of the need in GRASS to use a projected or geodetic system for any analysis. What is clear is that the outputs can be very different. For instance, notice that a degree in latitude does not have the same "size" as a degree in longitude.

Using equirectangular projection by default is probably hiding the distortion of that projection. This is not a conformal or equal area projection. IMHO, the user should select the projection if a cartesian system is needed.

If you are "converting" automatically EPSG:4326 into EPSG:4087 (https://epsg.org/crs_4087/WGS-84-World-Equidistant-Cylindrical.html) it should me mentioned somehow. But, what are you doing with other geographic systems like ETRS98 (EPSG:4258) or NAD83 (EPSG:4269)?

[florisvdh]

I think I initiated this subject at the GRASS GIS community meeting beginning of June, at the dinner table. Later I referred to EPSG:4087 in the grassgis/sprint gitter, as the CRS probably being used in display of EPSG:4326 data. Thanks @ldesousa for picking up and for your responses in gitter!

However I doubt that the default projection (world map case) should be set as EPSG:4087, since this would mean that the actual data are stored using this CRS. The location CRS (g.proj -p) must match the CRS of the data stored in that location IIUC.

projinfo EPSG:4087 -o WKT2_2019

``` $ projinfo EPSG:4087 -o WKT2_2019 WKT2:2019 string: PROJCRS["WGS 84 / World Equidistant Cylindrical", BASEGEOGCRS["WGS 84", ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)"], MEMBER["World Geodetic System 1984 (G730)"], MEMBER["World Geodetic System 1984 (G873)"], MEMBER["World Geodetic System 1984 (G1150)"], MEMBER["World Geodetic System 1984 (G1674)"], MEMBER["World Geodetic System 1984 (G1762)"], MEMBER["World Geodetic System 1984 (G2139)"], ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["World Equidistant Cylindrical", METHOD["Equidistant Cylindrical", ID["EPSG",1028]], PARAMETER["Latitude of 1st standard parallel",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8823]], PARAMETER["Longitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["False easting",0, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",0, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting (X)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (Y)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["Graticule coordinates expressed in simple Cartesian form."], AREA["World."], BBOX[-90,-180,90,180]], ID["EPSG",4087]] ```

My guess is that data in geographic CRSs are just displayed in the GRASS GUI using the 'World Equidistant Cylindrical' conversion method (equidistant cylindrical projection = equirectangular projection). So this would also hold for the other examples that @jjimenezshaw gives. But the data themselves will still be stored as geographic, which the GRASS commands should either be handling OK or else the user should indeed apply a projection (effectively creating a new location) before further processing.

At least, there should not be a problem with a location being geographic, hence with storing data using a geographic CRS in the GRASS database. I would have no problem with that being/remaining a default, although that's another facet of the discussion.

My original point was only the potential misunderstanding a user may have, since map rendering of geographic data always involves a projection while GRASS says that the CRS is (e.g.) EPSG:4326. So the user is looking at a projection while the data are geographic. Reminding the user of that fact somewhere, seems an appropriate way to deal with this IMHO. For example: "Data are geographic (EPSG:4326), rendered on the fly using the equirectangular projection method (EPSG:4087)."

[jjimenezshaw]

Minor detail.

"Data are geographic (EPSG:4326), rendered on the fly using the equirectangular projection method (EPSG:4087)."

should be "Data are geographic (EPSG:4326), rendered on the fly using the equirectangular projection method."

There are things like "methods" in EPSG, but EPSG:4087 is a CRS. In addition it would be the same string for every geographic system.

[florisvdh]

Oh yes, of course! Maybe indeed it's not needed to mention the CRS that results from the projection, although it's more informative to add which one, which stresses the CRS distinction. An alternative:

"Data are geographic (EPSG:4326), rendered using the equirectangular projection method (the map is in EPSG:4087)."

The used EPSG method:

METHOD["Equidistant Cylindrical",
            ID["EPSG",1028]]

[ldesousa]

Thank you @jjimenezshaw and @florisvdh for your insights. I would like to stress here that (most?) GRASS modules only apply to the Cartesian plane. You may well store your data in a geodetic CRS, but once you apply a module the Equirectangular projection is implicitly applied, say, to compute distances, vision basins, watersheds, etc. Map display is just one of such occasions.

Software will always benefit from being transparent, moreover considering the educational role GRASS intends to have. There is also the good example of PostGIS, where geodetic CRSs have their own data type and processing functions.

[florisvdh]

Thanks for adding these points @ldesousa. IMO, similarly to the proposed message regarding display, the modules should give a warning about this. Using equirectangular projection for planar calculations is often not good practice; for generic map display it may be an acceptable projection but often not so for calculations due to the distortions.

For example in R, the sf package (for vector data representation & processing) will by default use the S2 library in case of a geodetic CRS and no projections are performed (using GEOS only in case of projected CRS). And it will warn when doing planar calculations on data in a geographic CRS (assuming equirectangular projection):

> library(sf)
Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
> data_nad27 <- system.file("shape/nc.shp", package = "sf") |> read_sf()
> st_crs(data_nad27)$IsGeographic
[1] TRUE
> 
> sf_use_s2()
[1] TRUE
> data_buffered <- st_buffer(data_nad27, dist = 10)
> 
> sf_use_s2(FALSE)
Spherical geometry (s2) switched off
> data_buffered <- st_buffer(data_nad27, dist = 10)
dist is assumed to be in decimal degrees (arc_degrees).
Warning message:
In st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = endCapStyle,  :
  st_buffer does not correctly buffer longitude/latitude data

[wenzeslaus]

Many tools in GRASS GIS do indeed handle geographic CRS or give warning or error if they are not build for that. Some are too generic to make such decisions, but are doing their best to support user to do things right. I'm think r.mapcalc here, which does consider geographic CRS when calculating area of current cell in square meters with the area() function.

Code which handles geographic versus projected CRS is occurs throughout the code base. Here is a list of regular expressions I could think from top of my head and number their occurrences:

regular expression	number of occurrences
G_projection.*PROJECTION_LL	90
case PROJECTION_LL	8
G_begin_distance_calculations	51
Vect_line_geodesic_length	20
.*_geodesic_distance.*	52
G_begin_.*_area_calculations	21

What you wish for is already reality in GRASS GIS. Nevertheless, there might be some improvements needed and Grassy the Hungry Cow is always hungry for more code, but now, she just says:

  _________________________
/ I prefer loooooooongitude \
\            over latitude. /
  -------------------------
        \   ^__^
         \  (oo)\_______
            (__)\       )\/\
                ||----w |
_\|/__\|/__\|/_ ||     ||  _\|/__\|/_

[florisvdh]

Thanks for clarifying and illustrating this @wenzeslaus :+1:.

Then, the point may shrink to the original one, i.e. making user aware that map display is not EPSG:4326 (but EPSG:4087) even though the data are EPSG:4326. Or more generally, that 'geographic data are being displayed using the equirectangular projection'.