[R-sig-Geo] readOGR and proj4 string

[Agustin Lobo]

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

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

-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster