1

I'm trying to verify that my raster has the correct coordinate reference system represented in its WKT, as I had to manually re-assign the datum after importing it (due to the PROJ6 migration). However, the coordinate system identified in the WKT is EPSG 9122, which seems to be a transformation from NAD83(CRS)v2 to NAD83(CRS)v4 specific to Southwest Canada (https://epsg.io/9122) and not available on spatialreference.org. However, my raster covers Southwest Oregon, so I'm not sure where this came from.

Here is the spatial reference information from the metadata, which shows that the datum is NAD83 and only shows the horizontal coordinate system as latitude/longitude:

Spatial_Reference_Information:
  Horizontal_Coordinate_System_Definition:
    Geographic:
      Latitude_Resolution: 0.00001
      Longitude_Resolution: 0.00001
      Geographic_Coordinate_Units: Decimal degrees
    Geodetic_Model:
      Horizontal_Datum_Name: North American Datum of 1983
      Ellipsoid_Name: Geodetic Reference System 80
      Semi-major_Axis: 6378137.000000
      Denominator_of_Flattening_Ratio: 298.2572221
  Vertical_Coordinate_System_Definition:
    Altitude_System_Definition:
      Altitude_Datum_Name: North American Vertical Datum of 1988
      Altitude_Resolution: 0.001
      Altitude_Distance_Units: meters
      Altitude_Encoding_Method: Implicit coordinate

Data (including metadata): https://drive.google.com/drive/folders/1IhULan5vFPYZLzw6cDrlgklfYBGGhS-Z?usp=sharing

Here is the code:

> dir = "C:/Users/sarab/OneDrive/Documents/Thesis/R"
> setwd(dir)
> require(raster)
Loading required package: raster
Loading required package: sp
> require(sf)
Loading required package: sf
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
> #import elevation rasters
> elev <- raster("n43w124_13.img")
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
> #look at CRS as read into R
> st_crs(elev)
Coordinate Reference System:
  User input: +proj=longlat +ellps=GRS80 +no_defs 
  wkt:
GEOGCRS["unknown",
    DATUM["Unknown based on GRS80 ellipsoid",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1],
            ID["EPSG",7019]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
> #assign NAD83 as datum (per metadata) and look at CRS again
> crs(elev) <- "EPSG:4269"
> st_crs(elev)
Coordinate Reference System:
  User input: NAD83 (with axis order normalized for visualization) 
  wkt:
GEOGCRS["NAD83 (with axis order normalized for visualization)",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6269]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]] 


Did I assign the datum appropriately? Is there any way to find out more about "EPSG:9122" and why it has been assigned to this raster? It also seems that the ID for the datum has changed, from 4269 to 6269, so perhaps I'm missing something about the EPSG ID system.

1
  • If you want this looked at further, could you create a cropped or otherwise resized version of the raster with the same internal metadata, since a 499Mb file is going to take some time to download...
    – Spacedman
    Jan 2 at 12:48
2

EPSG codes re-use the same number for different purposes. Not only is 9122 the code for that transformation, it also the code for "degrees". From epsg.org:

enter image description here

Looking at the WKT, you can see that it is applied to the units of lat-long:

    ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]

So it is nothing to do with the EPSG:9122 transformation at all.

I've no idea why the EPSG use the same numbers for multiple types of things, but its usually clear from the context of the WKT which one is appropriate - but it can get obscured by the indenting and the plethora of square brackets.

Your confusion about 4269 and 6269 is the inverse problem - these two codes point to two different types of thing. epsg:4269 is a CRS (coordinate reference system) (and also an Extent, but you aren't talking about that) and 6269 is a Datum (and a Colombian CRS, but you aren't talking about that). The EPSG:4269 CRS is based on the ESPG:6269 Datum. The Datum defines things like the ellipsoid shape and the 0 lat and long points, a CRS puts other things on the Datum such as defining which way is up, or what units of measurement to use, or whether the CRS is projected or is lat-long.

5
  • epsg::number are just short hand, the full identifiers identify the type of thing
    – nmtoken
    Jan 2 at 16:19
  • well yes, "epsg:1234" is clearly a shorthand for something, being meaningless in itself. Why the EPSG use the same shorthand for two different things is a bit odd - its not like they are going to run out of unique numbers...
    – Spacedman
    Jan 2 at 17:09
  • Well the GML/WKT gives the context for the abridged version, otherwise use the full identifier
    – nmtoken
    Jan 2 at 17:12
  • What's a "full identifier"? How does it discriminate between the two EPSG:9122 items? Is it the id in the GML? Like gml:id="epsg-uom-9122" and epsg-op-9122? I've not seen those used in text before.
    – Spacedman
    Jan 2 at 17:17
  • They used to use URNs ~ urn:ogc:def:objectType:authority:version:code, where objectType identified the type of thing, like crs, datum, axis, unit...
    – nmtoken
    Jan 2 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.