On Mon, 11 Oct 2010, Agustin Lobo wrote:
Casper,
Thanks, this is very useful, although does not solve my problem.
The list of datums known to proj4 is given by projInfo("datum") - untried,
please try yourself. If they are not given, it is because they are
ambiguous, so proj4 sensibly errs on the side of caution. You have to give
the towgs84= values yourself for your area.
Roger
What I would need is to define:
proj4string(lakes) <- CRS("+proj=utm +zone=35 +datum=ED87
+towgs84=-82.981,-99.719,-110.709,-0.5076,0.1503,0.3898")
according to
http://www.asprs.org/resources/grids/10-99-norway.pdf
My question is why +datum=WGS84 works but +datum=ED87 or +datum=ED50
do not work?
Using A <- make_EPSG() I've been able to find:
A$note[grep("ED",A$note)]
A[160,]
? code ? note ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? prj4
160 4231 # ED87 +proj=longlat +ellps=intl +no_defs
but this is a lonlat projection, what I need is a UTM projection.
I've also found:
? ? code ? ? ? ? ? ? ? ? ?note
? ? prj4
2540 23035 # ED50 / UTM zone 35N +proj=utm +zone=35 +ellps=intl
+units=m +no_defs
and have used:
proj4string(lakes) <-CRS(paste("+init=epsg:23035",
"+towgs84=-82.981,-99.719,-110.709,-0.5076,0.1503,0.3898"))
lakesll <- spTransform(lakes, CRS("+proj=longlat +datum=WGS84"))
writeOGR(lakesll,dsn="Kolalakesll",layer="Kolalakesll",driver="ESRI
Shapefile")
but lakesll is shifted if overlaid on an image with lonlat and WGS84
"projection" and
I suspect the reason is using ED50 instead of ED87.
Note1. lakes comes from:
library(StatDA)
library(maptools)
xy <- kola.background$lakes
names(xy) <- c("x", "y")
lakes <- map2SpatialLines(xy)
lakesdata <- data.frame(ID=1:131,lakename=rep("",131),stringsAsFactors=F)
row.names(lakesdata) <- 1:131
lakes <- SpatialLinesDataFrame(sl=lakes, data=lakesdata, match.ID = TRUE)
Note2: sessionInfo()
R version 2.11.1 (2010-05-31)
i486-pc-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 ? ? ? ? ?LC_NUMERIC=C
LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 ? ? ? ?LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 ? ? ? ? ?LC_NAME=en_US.UTF-8
LC_ADDRESS=en_US.UTF-8
[10] LC_TELEPHONE=en_US.UTF-8 ? ? ?LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=en_US.UTF-8
attached base packages:
[1] stats ? ? graphics ?grDevices utils ? ? datasets ?methods ? base
other attached packages:
[1] rgdal_0.6-27 raster_1.1.2 sp_0.9-62 ? ?rkward_0.5.4
loaded via a namespace (and not attached):
[1] grid_2.11.1 ? ? lattice_0.19-11 tools_2.11.1
Thanks!
Agus
2010/10/11 caspar hallmann <caspar.hallmann at gmail.com>:
HI Agus,
see
A=make_EPSG()
A
in rgdal
Caspar
On Sun, Oct 10, 2010 at 10:36 PM, Agustin Lobo <alobolistas at gmail.com>
wrote:
CRS("+proj=utm +zone=32 +datum=WGS84")
CRS arguments:
?+proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
is correct, the same with ED50 is not:
CRS("+proj=utm +zone=32 +datum=ED50")
Error in CRS("+proj=utm +zone=32 +datum=ED50") :
?unknown elliptical parameter name
Is there any way to get the list of available datum acronyms?
Thanks
Agus