Hi Dominik,
If you use the gdalUtils package there is no significant difference in
the results using CLI or R:
library(gdalUtils)
gdaltransform(s_srs="+proj=longlat +datum=WGS84", t_srs="+proj=lcc
+lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0 +y_0=0
+ellps=sphere +a=6370000 +b=6370000 +units=m
+no_defs",coords=matrix(c(-112.25,33.0), ncol = 2))
[,1] [,2] [,3]
[1,] -1306676 -629522.5 0
If I use spTransform I can reproduce your results:
loc <- c(-112.25,33.0)
loc <- data.frame(matrix(unlist(loc), nrow=1,
byrow=T),stringsAsFactors=FALSE)
colnames(loc)<-c("lon","lat")
coordinates(loc)<- ~lon+lat
proj4string(loc)<- CRS("+proj=longlat +datum=WGS84 +no_defs")
loc<-spTransform(loc,CRS('+proj=lcc +lat_1=28 +lat_2=50
+lat_0=39.70001220694445 +lon_0=-98 +x_0=0 +y_0=0 +ellps=sphere
+a=6370000 +b=6370000 +units=m +no_defs'))
loc@coords
lon lat
1 -1306471 -629424
I suggest to focus on the sptransform() function
cheers Chris
> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.3 LTS
locale:
[1] LC_CTYPE=de_DE.UTF-8 LC_NUMERIC=C LC_TIME=de_DE.UTF-8
LC_COLLATE=de_DE.UTF-8
[5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=de_DE.UTF-8
LC_PAPER=de_DE.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] gdalUtils_2.0.1.7 raster_2.5-2 sp_1.2-2
loaded via a namespace (and not attached):
[1] rgdal_1.1-3 tools_3.2.3 Rcpp_0.12.3 R.methodsS3_1.7.0
codetools_0.2-14 grid_3.2.3
[7] iterators_1.0.8 foreach_1.4.3 R.utils_2.2.0
R.oo_1.19.0 lattice_0.20-33
Am 15.02.2016 um 21:37 schrieb Dominik Schneider:
Hi,
I'm struggling to use a custom projection. I am seeing differences with
someone using python proj4 bindings and when I compared my R results with
my commandline results I got even more confused. the coordinate
transformation is different for the two different methods.
could someone explain to me which one is wrong and why?
Thanks
R:
e=c(-112.25,-104.125,33,43.75)
box=as(extent(e),'SpatialPolygons')
proj4string(box)='+proj=longlat +datum=WGS84'
pstring='+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98
+x_0=0 +y_0=0 +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs'
xmin(spTransform(box,CRS(pstring)))
## [1] -1306471
ymin(spTransform(box,CRS(pstring)))
## [1] -713442.3
commandline:
iMac:~ $ gdaltransform -s_srs "+proj=longlat +datum=WGS84" -t_srs
"+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0
+y_0
=0 +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs"
-112.25 33.
-1306675.75472246 -629522.472322824 0
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
