read prj file to CRS

Roger Bivand a ?crit :

On Thu, 17 May 2007, Patrick Giraudoux wrote:

  

Dear all,

Definitely, I am not that clever with projection systems: I would like 
to read a *.prj file (part of a set of ESRI shapefiles) and get it in a 
CRS compatible format (handled by the  PROJ.4 projection system) and 
possibly added to sp objects where needed. What would would be the 
simplest way to do that? I am messing with rgdal and sp without result...

    

Yes, good idea, given that the reverse is provided in showWKT() in rgdal. 
As things are, you'll need to read it as part of a shapefile - readOGR() 
will read it when given a shapefile to read, or readGDAL() with an Arc 
Ascii grid, driver AAIGrid. I'll look at putting a simple *.prj reader 
into rgdal, but for now giving an arbitrary shapefile or Ascii grid the 
same name should do it.

Hope this helps,

Roger

  

OK, thanks. I went through it with readOGR(). I was messing with the 
readOGR documentation not understanding about how to read a shapefile 
specifically. Googling on the internet I found an example 
(https://stat.ethz.ch/pipermail/r-sig-geo/2007-February/001761.html) 
which made me make:

mymap<-readOGR("Mth030607.shp", layer="Mth030607")

which worked perfectly.

The trick was to understand that the first argurment 'dsn' expected a 
full shapefile name (before I tried the name "ESRI Shapefile" given by 
ogrDrivers(), and other deadlocks), and the second, 'layer', the file 
name truncated without its suffix. May I suggest to write an example 
about shapefiles in the readOGR() doc ?

Now there is no trouble to extract de CRS

getSlots(class(mymap))
        data   coords.nrs       coords         bbox  proj4string
"data.frame"    "numeric"     "matrix"     "matrix"        "CRS"

mymap at proj4string
CRS arguments:
 +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666667 
+k_0=0.99987742 +x_0=600000
+y_0=2200000 +a=6378249.2 +b=6356514.999904194 +units=m +no_defs

or more simply:

 > proj4string(mymap)

[1] " +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666667 
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 
+b=6356514.999904194 +units=m +no_defs"

Actually, I was used to WRITE *.prj files like that:

myproj<-showWKT(proj4string(mymap))
cat(myproj, file="Mth030607.prj")

Now I can close the circle writing AND reading... I am quite lazy with 
projection files (writing is source of many mistakes), so it will be 
easy to pick them up from already existing *.prj files and to create and 
add CRS to spatial objects and write the corresponding *.prj when 
shapefiles are written (eg with readShapePoly(), etc.). No longer need 
the ESRI ArcToolbox.

Thanks,

Patrick