[R-sig-Geo] Bearing angle of UTM projected SpatialLines
Roger Bivand
Roger.Bivand at nhh.no
Thu May 26 12:28:04 CEST 2016
On Wed, 25 May 2016, Eduardo Diez wrote:
> I forgot to add the complete coordinates for the line:
>
>> l1 at lines[[1]]@Lines[[1]]@coords
> [,1] [,2]
> [1,] 642236.2 6197243
> [2,] 643005.6 6197948
>
>> l1.geo at lines[[1]]@Lines[[1]]@coords
> [,1] [,2]
> [1,] -61.45336 -34.35639
> [2,] -61.44511 -34.34994
#This was:
library(sp)
l1c <- cbind(c(642236.2, 643005.6), c(6197243, 6197948))
lgc <- cbind(c(-61.45336, -61.44511), c(-34.35639, -34.34994))
l1 <- SpatialLines(list(Lines(list(Line(l1c)), "1")),
proj4string=CRS("+proj=utm +zone=20 +south +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0"))
l1.geo <- SpatialLines(list(Lines(list(Line(lgc)), "1")),
proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"))
#but the lack of precision (I asked for binary) in the coordinates means
#that Arc's:
#Angle given by Linear Directional Mean for l1:
#47.53076
#becomes:
ds <- apply(coordinates(l1)[[1]][[1]], 2, diff)
atan2(ds[1], ds[2])*180/pi
atan2(ds[1], ds[2]-0.7343)*180/pi
#because the rounded difference in Northings is off by 0.7m.
#The spherical/ellipsoidal calculations are:
library(maptools)
gzAzimuth(from = l1.geo at lines[[1]]@Lines[[1]]@coords[1,], to =
l1.geo at lines[[1]]@Lines[[1]]@coords[2,], type="snyder_sphere")
library(geosphere)
bearing(p1 = l1.geo at lines[[1]]@Lines[[1]]@coords[1,], p2 =
l1.geo at lines[[1]]@Lines[[1]]@coords[2,], sphere=TRUE)
#equal as spherical, with this differeing a little depending where you are
#on the globe:
bearing(p1 = l1.geo at lines[[1]]@Lines[[1]]@coords[1,], p2 =
l1.geo at lines[[1]]@Lines[[1]]@coords[2,], sphere=FALSE)
#obviously one is spherical, the other ellipsoidal. Neither are planar, so
#are not the same thing anyway.
If anyone feels like contributing a function to accept objects rather than
raw coordinates for gzAzimuth, and restricting to geographical
coordinates, adding a planar direction function, and an LDM function for
Lines and SpatialLines objects, or any of these, I'd be grateful.
Hope this clarifies,
Roger
>
> Thanks
>
> 2016-05-25 16:09 GMT-03:00 Eduardo Diez <eduardodiez at gmx.com>:
>
>> Here are some details and example:
>>
>>> sessionInfo()
>> R version 3.2.5 (2016-04-14)
>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>> Running under: Windows 7 x64 (build 7601) Service Pack 1
>>
>> locale:
>> [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
>> States.1252
>> [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
>>
>> [5] LC_TIME=English_United States.1252
>>
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> other attached packages:
>> [1] RPyGeo_0.9-3 RSAGA_0.94-5 plyr_1.8.3
>> gstat_1.1-3
>> [5] shapefiles_0.7 foreign_0.8-66 rgdal_1.1-8
>> maptools_0.8-39
>> [9] geosphere_1.5-1 ggplot2_2.1.0 raster_2.5-2
>> rgeos_0.3-19
>> [13] sp_1.2-3 RevoUtilsMath_3.2.5
>>
>> loaded via a namespace (and not attached):
>> [1] Rcpp_0.12.4 maps_3.1.0 munsell_0.4.3 colorspace_1.2-6
>> lattice_0.20-33
>> [6] FNN_1.1 xts_0.9-7 tools_3.2.5 parallel_3.2.5
>> grid_3.2.5
>> [11] gtable_0.2.0 intervals_0.15.1 digest_0.6.9 mapproj_1.2-4
>> labeling_0.3
>> [16] scales_0.4.0 spacetime_1.1-5 zoo_1.7-12
>>
>> Spatial line: l1
>>> str(l1)
>> Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
>> ..@ data :'data.frame': 1 obs. of 2 variables:
>> .. ..$ OBJECTID : int 1
>> .. ..$ SHAPE_Leng: num 1043
>> .. ..- attr(*, "data_types")= chr [1:2] "N" "F"
>> ..@ lines :List of 1
>> .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
>> .. .. .. ..@ Lines:List of 1
>> .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
>> .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 642236 643006 6197243
>> 6197948
>> .. .. .. ..@ ID : chr "0"
>> ..@ bbox : num [1:2, 1:2] 642236 6197243 643006 6197948
>> .. ..- attr(*, "dimnames")=List of 2
>> .. .. ..$ : chr [1:2] "x" "y"
>> .. .. ..$ : chr [1:2] "min" "max"
>> ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
>> .. .. ..@ projargs: chr "+proj=utm +zone=20 +south +datum=WGS84 +units=m
>> +no_defs +ellps=WGS84 +towgs84=0,0,0"
>>
>> Spatial line in geographic coordinates: l1.geo
>>> str(l1.geo)
>> Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
>> ..@ data :'data.frame': 1 obs. of 2 variables:
>> .. ..$ OBJECTID : int 1
>> .. ..$ SHAPE_Leng: num 1043
>> .. ..- attr(*, "data_types")= chr [1:2] "N" "F"
>> ..@ lines :List of 1
>> .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
>> .. .. .. ..@ Lines:List of 1
>> .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
>> .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] -61.5 -61.4 -34.4 -34.3
>> .. .. .. ..@ ID : chr "0"
>> ..@ bbox : num [1:2, 1:2] -61.5 -34.4 -61.4 -34.3
>> .. ..- attr(*, "dimnames")=List of 2
>> .. .. ..$ : chr [1:2] "x" "y"
>> .. .. ..$ : chr [1:2] "min" "max"
>> ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
>> .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs
>> +ellps=WGS84 +towgs84=0,0,0"
>>
>> Run of maptools::gzAzimuth:
>>> gzAzimuth(from = l1.geo at lines[[1]]@Lines[[1]]@coords[1,],
>> + to = l1.geo at lines[[1]]@Lines[[1]]@coords[2,])
>> [1] 46.52677
>>
>> Run of geosphere::bearing:
>> bearing(p1 = l1.geo at lines[[1]]@Lines[[1]]@coords[1,],
>> + p2 = l1.geo at lines[[1]]@Lines[[1]]@coords[2,])
>> [1] 46.65786
>>
>> Angle given by Linear Directional Mean for l1:
>> 47.53076
>>
>> For l1.geo (in Lat/Lon) the ArcGIS' tool gives a warning that it it should
>> be projected but calculates this angle:
>> WARNING 000916: The input feature class does not appear to contain
>> projected data.
>> 51.94714
>>
>> I can provide further information or files but i think this can be enough
>> for now.
>> Thanks,
>> Regards
>>
>> 2016-05-25 4:32 GMT-03:00 Roger Bivand <Roger.Bivand at nhh.no>:
>>
>>> On Wed, 25 May 2016, Eduardo Diez wrote:
>>>
>>> Dear everyone,
>>>> I'm used to calculating the compass angles (clockwise from due North) of
>>>> line features projected in UTM using the tool Linear Directional Mean
>>>> <
>>>> http://resources.esri.com/help/9.3/arcgisengine/java/gp_toolref/spatial_statistics_tools/how_linear_directional_mean_spatial_statistics_works.htm
>>>>>
>>>> from ArcGIS.
>>>> I could find functions for performing a similar task but taking Origin ->
>>>> Destination points and only in Lat/Lon (geographic coordinates) namely:
>>>> geosphere::bearing and maptools::gzAzimuth. (Somehow they give different
>>>> results for the same set of points: around 0.05 degrees, maybe because of
>>>> the trigonometry)
>>>> Because of the nature of my work it has to be in UTM.
>>>>
>>>
>>> Interesting question. Could you please provide a test set of pairs of
>>> coordinates with the fully qualified UTM projection (which datum in
>>> particular, and if WGS84, which version - there are several), the output
>>> from ArcGIS (which version), and the geographical coordinates you see in
>>> ArcGIS (to remove the possibility that it is a projection issue on the R
>>> side). They don't need to mean anything, only that they can show us what
>>> the problem is in code and numeric terms. A link to an R file and an RDS
>>> file would be helpful.
>>>
>>> Roger
>>>
>>> Perhaps someone skillfull in trigonometry can figure it out but i'm not
>>>> the
>>>> case.
>>>> Is there a way to do this in R?
>>>>
>>>> Thanks,
>>>> Eduardo
>>>>
>>>> [[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
>>>>
>>>>
>>> --
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>> e-mail: Roger.Bivand at nhh.no
>>> http://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>> http://depsy.org/person/434412
>>>
>>
>>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412
More information about the R-sig-Geo
mailing list