Bearing angle of UTM projected SpatialLines

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

Bearing angle of UTM projected SpatialLines

Thread (5 messages)