readOGR and proj4 string

On Thu, 10 Apr 2008, Agustin Lobo wrote:

Thanks,

This is what I'm doing:

 absUTMpolysHABS2 <- readOGR("C:/ALOBO/Lidia",layer="test_TNT")
 t1 <- CRS(paste("+proj=utm +zone=31 

+ellps=intl","+towgs84=-87.0,-98.0,-121.0,0.0,0.0,0.0,0.0"))

 proj4string(absUTMpolysHABS2) <-t1
 proj4string(x) <-t1

As my goal is

 a <- (joinPolys(SpatialPolygons2PolySet(x), 

SpatialPolygons2PolySet(absUTMpolysHABS2)))

the critical issue is that both objects have the
same proj4string. But I'm saving absUTMpolysHABS2 again
as shp and will look at it carefully to clarify
the problem of the ED50 definition.

But I'm getting this problem, which seems to be related to the
proj4string:

 a1 <- SpatialPolygons2PolySet(absUTMpolysHABS2)
 a2 <- PolySet2SpatialPolygons(a1)

Error in CRS(p4s) : invalid UTM zone number

Please ensure that t1 looks sane, spaces between tag=value pairs.

PBS have changed the way they store the UTM zone - it is now a top level 
attribute, but was an attribute of an attribute:

t1 <- CRS(paste("+proj=utm +zone=31 +ellps=intl",
   "+towgs84=-87.0,-98.0,-121.0,0.0,0.0,0.0,0.0"))
proj4string(x) <- t1
a1 <- SpatialPolygons2PolySet(x)
attr(a1, "projection")
attr(attr(a1, "projection"), "zone") <- attr(a1, "zone")
attr(a1, "projection")
a2 <- PolySet2SpatialPolygons(a1)

I'll update the function in due course.

Hope this helps,

Roger

 summary(a1)

PolySet

Records         : 95995
  Contours      : 1595
    Holes       : 0
  Events        : NA
    On boundary : NA

Ranges
  X             : [412000, 466000]
  Y             : [4584000, 4624000]

Attributes
  Projection    : UTM
  Zone          : 31

Extra columns   :

Agus

Roger Bivand escribi?:

 On Thu, 10 Apr 2008, Roger Bivand wrote:

 On Thu, 10 Apr 2008, Agustin Lobo wrote:

 The manual page of readOGR states:

 p4s     PROJ4 string defining CRS, if default NULL, the 

value is read > >  from

 the OGR data set

 then, if a  *.prj file is present for a *.shp,

 why is the proj4string of the resulting SpatialPolygonsDataframe
 set to NA? is this a general behaviour or am I doing someting
 wrong?
 Puzzling. Could you make your test file available for me to check?

 OK, thanks. I can reproduce the problem. PROJ.4 does not support
 projection "names" as such, they are often not unique (same name but
 different parameters. In addition, ED50 is a collection of 
standards, not
 one standard. I suggest entering:

 t1 <- CRS(paste("+proj=utm +zone=31 +ellps=intl",
  "+towgs84=-87.0,-98.0,-121.0,0.0,0.0,0.0,0.0")

 manually. http://www.asprs.org/resources/grids/, July 200, gives 
different
 +towgs84= values. The ones in the file are prepended by a comma, which
 looks odd, but removing it doesn't help in parsing the file.

 CRS("+init=epsg:23031")

 is correct, but as usual with EPSG, no +towgs84= is given for ED50,
 because they vary so much.

 So the underlying code in readOGR is erring on the side of caution - if
 the *.prj file cannot be interpreted unequivocally, return NA.

 I set the value above:

 x <- readOGR(".", "test_TNT")
 proj4string(x) <- t1
 x1 <- spTransform(x, CRS("+proj=longlat +datum=WGS84"))
 writeOGR(x1, "x1.kml", "x1", driver="KML")

 and viewed in Google Earth, things looked more-or-less OK. You could 
try
 the +towgs84 parameters from Grids & Datums and see if they fit 
better on
 the ground.

 Hope this helps,

 Roger

 Roger

 absUTMpolysHABS2 <- 

readOGR("C:/ALOBO/Lidia",layer="test_TNT")

 where I have:
 test_TNT.avl

 test_TNT.dbf
 test_TNT.prj
 test_TNT.shp
 test_TNT.shx

 with test_TNT.prj:
 

PROJCS["ED50_/_UTM_zone_31N_(CM_3E)",GEOGCS["ED50_/_Geographic",DATUM["D_European_1950",SPHEROID["International_1924",6378388.0,297.0]TOWGS84[,-87.0,-98.0,-121.0,0.0,0.0,0.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Decimal_Degree",0.01745329251994330]],PROJECTION["Transverse_Mercator"],PARAMETER["Latitude_Of_Center",0.0],PARAMETER["Longitude_Of_Origin",3.0],PARAMETER["Scale_Factor",0.9996000000],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],UNIT["Meter",1.0]] 

 I get:

 class(delme)

 [1] "SpatialPolygonsDataFrame"
 attr(,"package")
 [1] "sp"

 str(delme,max.level=2)

 Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
   ..@ data       :'data.frame': 1555 obs. of  16 variables:
   ..@ polygons   :List of 1555
   ..@ plotOrder  : int [1:1555] 148 143 792 740 209 335 895 619 

1127

 1160 ...
   ..@ bbox       : num [1:2, 1:2]  412000 4584000  466000 4624000
   .. ..- attr(*, "dimnames")=List of 2
   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots

 delme at proj4string

 CRS arguments: NA

 Agus