Skip to content

Bearing angle of UTM projected SpatialLines

5 messages · Eduardo Diez, Roger Bivand

Eduardo Diez
# 
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.
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
Roger Bivand
#
On Wed, 25 May 2016, Eduardo Diez wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Eduardo Diez
    

    

      
      #
    

  


  

      
Here are some details and example:

            

      
      
      
    

        
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

            

      
      
      
    

        
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

            

      
      
      
    

        
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:

            

      
      
      
    

        
+                   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>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Eduardo Diez
    

    

      
      #
    

  


  

      
I forgot to add the complete coordinates for the line:

            

      
      
      
    

        
[,1]    [,2]
[1,] 642236.2 6197243
[2,] 643005.6 6197948

            

      
      
      
    

        
[,1]      [,2]
[1,] -61.45336 -34.35639
[2,] -61.44511 -34.34994

Thanks

2016-05-25 16:09 GMT-03:00 Eduardo Diez <eduardodiez at gmx.com>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 25 May 2016, Eduardo Diez wrote:

        

            

      
      
      
    

        
#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