2

I am working with the North America Forest Dynamics dataset. As mentioned it uses the USA Albers Equal Area USGS version projection, which is primarily different from the USA Albers Equal Area projection in the measure of the latitude of origin (former is 23 and latter is 37.5).

The datum is clearly read by Arc as D_NorthAmerica 1983. which is the datum. However, this datum information is not read by R (using raster package). Hence, resulting analysis and written rasters have unknown datum when processing is done in R.

Interestingly, the famous spatial reference website does not have a Proj4 string info for this, which I could manually input into the crs() command. But on the same website, the .Prj file can be downloaded and on opening this ile as a text pad file a datum i.e. North America 1983 datum and hence transformation is mentioned.

Why is R not recognizing the datum around which the transformation has to happen for this projection?

1

I created a small TIFF from the data selector on that site. gdalinfo reports this:

PROJCS["unnamed",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

and R reads it and reports this:

> r = raster("./sdat_1290_1_20180804_101230826.tif")
> crs(r)
CRS arguments:
 +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

The "datum" info is encoded in the +ellps=GRS80 and the +towgs84=[zeroes] parameters. What's missing is the "name" in the DATUM section, namely, "North_American_Datum_1983", but I don't think these are official text strings in any sense. Some notes here on that: How are ESRI WKT projections different from OGC WKT projections?

So R is reading the correct data, and projections and transformations should work correctly.

  • 1
    As you mention, the crs of USGS Albers Equal Area Conic comes out to be exactly what you have written, specifically the +towgs84=0,0,0,0,0,0,0 which contains the datum info. Using spTransform, I reprojected a shp,in WGS 1984, to USGS projection. In R, I see that the new projected shapefile has the USGS albers projection string. Once I write it using writeOGR,import it back into R or look at it in QGIS, the +towgs84=0,0,0,0,0 is missing in the proj string and it comes up as D_unknown datum in ArcMap. Is there something wrong in the writeOGR()? I can edit my question above & include my code. – tg110 Aug 6 '18 at 16:17
  • curious if anyone has any updates since last year? I have come across this projection in a lot of NOAA hydrology data, and would like to be able to trade processed data back and forth between ArcMap and R. – dbo Jul 20 at 5:23

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.