Re: [R-sig-Geo] different projection transformation R and gdal commandline

[Dominik Schneider]

Oh, silly me. doing plot(spTransform(box,CRS(pstring))) and it's obvious
what's happening. The projection rotates the polygon such that for the
corners y1,x1 and y2,x2   y1 != y2. But ymin(spTransform(box,CRS(pstring)))
still gives you the smallest coordinate regardless of which corner it is.

Using spatial points its easier to see.

e=c(-112.25,-104.125,33,43.75)
box=as(extent(e),'SpatialPoints')
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=637 +b=637 +units=m +no_defs'
coordinates(spTransform(box,CRS(pstring)))
## x y
##[1,] -1306471.3 -629424.0
##[2,] -1122174.4  531025.6
##[3,]  -563451.2 -713442.3
##[4,]  -483968.2  458859.3




On Tue, Feb 16, 2016 at 9:27 AM, Dominik Schneider <
dominik.schnei...@colorado.edu> wrote:

> Hi Chris,
> Thanks for confirming this. I'm not surprised that gdalUtils gives the
> same answer as the gdal utilities - my understanding is that gdalUtils is
> basically the equivalent to calling the commandline utilities via system().
> I'm hoping that someone can shed light on spTransform since I use that a
> lot for transforming points and polygons.
> ds
>
>
> On Mon, Feb 15, 2016 at 2:50 PM, Chris Reudenbach <
> reudenb...@uni-marburg.de> wrote:
>
>> 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=637 +b=637 +units=m
>> +no_defs",coords=matrix(c(-112.25,33.0), ncol = 2))
>>  [,1]  [,2] [,3]
>> [1,] -1306676 -629522.50
>>
>>
>> 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=637
>> +b=637 +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-8LC_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=637 +b=637 +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=637 +b=637 +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
>>>
>>>
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] different projection transformation R and gdal commandline

[Dominik Schneider]

Hi Chris,
Thanks for confirming this. I'm not surprised that gdalUtils gives the same
answer as the gdal utilities - my understanding is that gdalUtils is
basically the equivalent to calling the commandline utilities via system().
I'm hoping that someone can shed light on spTransform since I use that a
lot for transforming points and polygons.
ds


On Mon, Feb 15, 2016 at 2:50 PM, Chris Reudenbach  wrote:

> 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=637 +b=637 +units=m
> +no_defs",coords=matrix(c(-112.25,33.0), ncol = 2))
>  [,1]  [,2] [,3]
> [1,] -1306676 -629522.50
>
>
> 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=637
> +b=637 +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-8LC_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=637 +b=637 +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=637 +b=637 +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
>>
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] different projection transformation R and gdal commandline

[Chris Reudenbach]

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=637 +b=637 +units=m 
+no_defs",coords=matrix(c(-112.25,33.0), ncol = 2))

 [,1]  [,2] [,3]
[1,] -1306676 -629522.50


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=637 +b=637 +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-8LC_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=637 +b=637 +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=637 +b=637 +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

Re: [R-sig-Geo] different projection transformation R and gdal commandline

[Dominik Schneider]

Sorry, forgot to add the details. rgdal and gdal 1.11.3 were installed from
kyngchaos.com

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-apple-darwin14.5.0 (64-bit)
Running under: OS X 10.11.3 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats graphics  grDevices utils datasets  methods
base

other attached packages:
 [1] rasterVis_0.37  latticeExtra_0.6-26 lattice_0.20-33
broom_0.4.0
 [5] tidyr_0.3.1 dplyr_0.4.3 doMC_1.3.4
 iterators_1.0.8
 [9] foreach_1.4.3   ncdf4_1.15  gridExtra_2.0.0
spdep_0.5-92
[13] Matrix_1.2-3ipred_0.9-5 MASS_7.3-45
RColorBrewer_1.1-2
[17] rgdal_1.0-7 stringr_1.0.0   ggplot2_2.0.0   plyr_1.8.3

[21] reshape2_1.4.1  raster_2.5-2sp_1.2-1
 gstat_1.1-0
[25] ProjectTemplate_0.6

loaded via a namespace (and not attached):
 [1] zoo_1.7-12   splines_3.2.3colorspace_1.2-6 spacetime_1.1-5
 [5] survival_2.38-3  prodlim_1.5.7hexbin_1.27.1DBI_0.3.1
 [9] lava_1.4.1   munsell_0.4.2gtable_0.1.2 codetools_0.2-14
[13] coda_0.18-1  psych_1.5.8  labeling_0.3 class_7.3-14
[17] xts_0.9-7Rcpp_0.12.2  scales_0.3.0 FNN_1.1
[21] deldir_0.1-9 mnormt_1.5-3 digest_0.6.8 stringi_1.0-1
[25] grid_3.2.3   tools_3.2.3  LearnBayes_2.15  magrittr_1.5
[29] lazyeval_0.1.10  assertthat_0.1   R6_2.1.1 rpart_4.1-10
[33] boot_1.3-17  intervals_0.15.1 nnet_7.3-11  nlme_3.1-122


Dominik Schneider
c 518.956.3978


On Mon, Feb 15, 2016 at 1:37 PM, Dominik Schneider <
dominik.schnei...@colorado.edu> wrote:

> 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=637 +b=637 +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=637 +b=637 +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] different projection transformation R and gdal commandline

[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=637 +b=637 +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=637 +b=637 +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