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@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





	[[alternative HTML version deleted]]


