[R-sig-Geo] readOGR and proj4 string
Agustin Lobo
Agustin.Lobo at ija.csic.es
Thu Apr 10 17:44:49 CEST 2008
It works fine now, thanks.
Then, back to the original problem of recovering
a SpatialPolygonsDataFrame from a PolySet
(which I get after joinPolys() ):
a1 <- SpatialPolygons2PolySet(x)
attr(a1, "projection")
attr(attr(a1, "projection"), "zone") <- attr(a1, "zone")
attr(a1, "projection")
a2 <- PolySet2SpatialPolygons(a1)
a3 <- SpatialPolygonsDataFrame(a2,x at data)
where a3 is identical to x, right?
Almost there. The problem with the ring
after joinPolys() remains, though. Also, once
I'll able to get the result of joinPolys() back to
SpatialPolygons I will have to subset and update
the original data.frame.
Agus
Roger Bivand escribió:
> 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
>>> > > > >
>>>
>>
>>
>
--
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
More information about the R-sig-Geo
mailing list