[R-sig-Geo] read prj file to CRS

[Roger Bivand]

On Thu, 17 May 2007, Patrick Giraudoux wrote:

> 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")

Or simpler:

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

where dsn= is a directory (data source name) containing one or more 
shapefiles. Because OGR supports many formats, it has to handle many 
different ways of defining the data sources, and that is why dsn= and 
layer= are split in this way.

writeOGR() uses similar mechanisms to readOGR(), and for most users should 
be prefered to the maptools shapefile functions. The exception is if you 
have difficulty in installing rgdal on your platform, because both 
Linux/Unix and Mac OSX need to link to external libraries.

Roger

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

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no