[R-sig-Geo] readOGR and proj4 string

[Roger Bivand]

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

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