5

I am working with a raster in the .e00 format, located here:

ftp://ftp.epa.gov/castnet/tdep/grids/n_dw/n_dw-2000.zip

I unzip, and then load this raster in R:

require(raster)
require(rgdal)
require(sp)
r <- raster('/path/to/n_dw-2000.e00')

Which works fine:

> r.d.2000
class       : RasterLayer 
dimensions  : 775, 1440, 1116000  (nrow, ncol, ncell)
resolution  : 4134.383, 4134.383  (x, y)
extent      : -2950369, 3003143, 115686.8, 3319834  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=0 +lat_2=29.5 +lat_0=45.5 +lon_0=0 +x_0=0 +y_0=-96 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : /Users/colin/Documents/Talbot Lab/FIA_2016_manuscript_EM-ndep/extracting_e00.files_dec2016/dry_dep/n_dw-2000.e00 
names       : n_dw.2000 
values      : 0.4092157, 72.94985  (min, max)

I then load some lat/long values I wish to extract the values from the raster for.

Here are 6 sites, with columns corresponding to longitude and latitude:

points <- structure(c(-86.030373, -85.990492, -85.847991, -87.954285, -87.674475, 
-85.138866, 46.414597, 46.137664, 42.258794, 46.567187, 45.730901, 
44.383864), .Dim = c(6L, 2L))
> points
          [,1]     [,2]
[1,] -86.03037 46.41460
[2,] -85.99049 46.13766
[3,] -85.84799 42.25879
[4,] -87.95428 46.56719
[5,] -87.67448 45.73090
[6,] -85.13887 44.38386

I need to project these lat/long values into the coordinate reference system (CRS) of my raster. To do this I first convert the lat/long to a SpatialPoints data frame:

points.2 <- SpatialPoints(points, proj4string = CRS("+init=epsg:4326"))

I then re-project these points into the CRS of the raster using spTransform

points.3 <- spTransform(points.2, projection(r))  

Finally, I would use the extract function to pull values from the raster for all of my points using:

out <- extract(r, points.3)

However this only returns NA's, to my dismay:

> out
     [,1]
[1,]   NA
[2,]   NA
[3,]   NA
[4,]   NA
[5,]   NA
[6,]   NA

I think the problem is how I am re-projecting the points, as if I plot the points vs. the raster they seem to be in a different orientation.

Furthermore, if I go to plot the points on otp of the raster I do not have any luck. I am unsure what specifically is going wrong however.

5

The coordinate system information is incorrect. You can usually find the information in a .e00 file. Here's what it says:

spheroid = GRS80
central meridian/longitude of origin = -96.0
standard parallel 1 = 29.5
standard parallel 2 = 45.5
latitude of origin = 23.0
false easting/false northing = 0.0

so the PROJ.4 string should be:

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23.0 +lon_0=-96.0 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

To change the CRS of your raster use the crs command:

crs(r) <- '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23.0 +lon_0=-96.0 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'

Then continue re-projecting the coordinates of your sites into the CRS of the raster with the code you have written.

Getting NaNs ("NA") back is sometimes a sign that the input coordinates can't be unprojected/reprojected which in turn can mean that the coordinate system information is incorrect. I've seen those parameter values before, but they weren't assigned to the parameters that I expected.

4
  • thanks so much. How did you obtain the correct CRS from the raster, and any hunch why it prints the 'wrong' CRS when loaded in R?
    – user234105
    Dec 29 '16 at 20:07
  • also- can you edit the code posted to account for this and extract correctly?
    – user234105
    Dec 29 '16 at 20:09
  • 1
    If an .e00 file has a defined coordinate system you can search for it in the file by looking for "PRJ". The description of the parameter follows the parameter value. No idea why it's being misinterpreted. It's not a mix-up like thinking the value follows the parameter description.
    – mkennedy
    Dec 29 '16 at 20:12
  • I'm not an R programmer. Maybe someone else will know how to fix it.
    – mkennedy
    Dec 29 '16 at 20:13

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.