[R-sig-Geo] readOGR and proj4 string

[Roger Bivand]

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