1

This is a follow-up post to a recent question on StackOverflow. I'm now facing a GIS problem rather than a coding one.

I'm working with a DEM from this source, which has the following projection information: +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0.

I need to projectRaster(mydem, worldclim) to fit worldclim data from this source, which has the following projection information: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0.

However, when I reproject mydem using worldclim as a template, I run into major alignment issues. mydem is a map of North America in Lambert azumithal equal-area and worldclim is the full globe in lat/lon, so I'm wondering if there is something going on extent-wise that might be causing the problem.

Here is an image of the original mydem: mydem

Here is an image of the original worldclim layer: worldclim

Here is an image of the reprojected DEM, calculated using

mydem.lonlat <- projectRaster(mydem, worldclim)

reprojected DEM

And here are the mydem.lonlat and worldclim layers overlain with one another: Overlay

I followed Dr. Hijman's approach in this post but I don't think rotation is the issue here. There are a few other related questions here and here but they are unanswered and still don't address my exact issue.

extent(mydem.lonlat) and extent(worldclim) are identical, as are res() and proj4string() for both layers... so they definitely align in plotting space; there's just clearly an issue with geographic space.

2

R read a weird projection. The original one is (from documentation attached in raster data folder):

PROJCS["NAD_1983_Lambert_Azimuthal_Equal_Area",GEOGCS 
      ["GCS_Sphere_ARC_INFO",DATUM["D_Sphere_ARC_INFO",SPHEROID 
      ["Sphere_ARC_INFO",6370997.0,0.0]],PRIMEM["Greenwich",0.0], 
      UNIT["Degree",0.0174532925199433]],PROJECTION 
      ["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",0.0], 
      PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian", 
      -100.0],PARAMETER["Latitude_Of_Origin",45.0],UNIT["Meter",1.0]]

So, this data projected US National Atlas Equal Area (EPSG:2163), not in WGS Datum as R read it. Just correct projection and will work:

library(raster)

mydem <- raster('~/Downloads/Elevation_GRID/NA_Elevation/data/NA_Elevation/na_elevation/hdr.adf')

projection(mydem) <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

worldclim <- getData('worldclim', var = 'tmin', res = 10)

mydem.lonlat <- projectRaster(mydem, worldclim)

plot(mydem.lonlat)

enter image description here

  • Terrific, thanks for your help. To follow up, why would R interpret the projection information incorrectly in the first place? It would make more sense to me if R just gave some sort of error rather than assume a particular set of projection parameters. Also, how did you go from the information in the documentation to a concise proj4string? Do you have any good relevant sources for translating between the original projcs to R's proj4string? – JepsonNomad Aug 1 '18 at 5:22
  • 1
    The projection strings are not standardized, which leads to many different variants. It should however lead to an error or that the program uses the constants and ellipsoid correctly, even when using another EPSG. But it seems that the information in the DEM was already wrong (through a wrong export maybe). Having such problems quite often i therefore tend to not use the projection information of the file but set the projection in the software itself based on the one in the documentation (if there is any). A good point for comparing the projection strings is for example: epsg.io/2163 – Matte Aug 1 '18 at 11:00
  • 1
    I just recognize projection. You can try to translate parameters to a proj.4 format with osr, but this projection cannot be translated automatically to a proj.4 format because some parameters are directly from ArcGIS, like "GCS_Sphere_ARC_INFO", "D_Sphere_ARC_INFO" or "Sphere_ARC_INFO". Because of this, R misread projection – aldo_tapia Aug 1 '18 at 13:44

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.