[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=6370000 +b=6370000 +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.schneider at 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 <
> reudenbach at 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=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 at 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 at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>
>

	[[alternative HTML version deleted]]