Defining datum for Lambert Azimuthal Equal Area and converting to geographic coordinates

Question

I posted this question to a specialised mailing list for spatial data analysis using R here but didn't recieve a response.

I am working with shapefiles and rasters in R, and I am unsure if I am specifying my projection information correctly (which I need to make sure I am as I want to reproject from one coordinate system to another).

I have a shapefile (from the Hydro1K data). The projection information in the readme is described as:

Projection used:  Lambert Azimuthal Equal Area
    Units = meters
    Pixel Size = 1000 meters
    Radius of Sphere of Influence = 6,370,997 meters
    Longitude of Origin = 20 00 00E
    Latitude of Origin = 5 00 00N
    False Easting = 0.0
    False Northing = 0.0

Looking at information in the readme and the associated .prj file I set the projection information using proj arguments like so:

proj <- CRS("+proj=laea +lon_0=20 +lat_0=5 +ellps=sphere") #the ellps 'sphere' has a radius of 6370997.0m
projection( myshpfile ) <- proj # In R the projection function allows you to set projection info for a file

However, I was wondering, do I also need to specify a datum for this projection? In the .prj file it says:

DATUM["D_Sphere_ARC_INFO",SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]]

But I am not sure what datum this refers to? Does anyone know? Or can tell me where to find out?

I also have some data which is in geographic coordinates. The .prj file of this data states:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Which I think is standard geographic WGS84 projection? Am I correct in thinking that the proj arguments for this projection should be:

projgeo <- CRS("+proj=latlong +datum=WGS84 +ellps=WGS84")
projection(mygeoshpfile) <- projgeo`

Then to convert mygeoshpfile to the same projection as myshpfile, in R I can use spTransform like this:

mygeoshpfile2 <- spTransform(mygeoshpfile, proj)

I want to make sure I am accurately converting from one projection to the other. Do I need to specify further proj parameters?

Thanks

Answer 1

This may be a silly question, but it looks like you have a .prj file, so why not just use rgdal to read it directly?

library (rgdal)

hydro<-readOGR(".","myshapefile") #read in file
proj4string(myshapefile) #get a printout of your file's projection in proj4 format

as for the second part of your question. Think about what it would mean to reproject a raster. Raster is just a grid. Reprojection would change the regularity on which the entire raster structure is based. Much better to convert going the other way using

myshape2<-spTransform(myshapefile,projgeo)