Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
add comment

1 Answer

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)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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