[R-sig-Geo] Bearing angle of UTM projected SpatialLines

Eduardo Diez eduardodiez at gmx.com
Wed May 25 21:09:37 CEST 2016 

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
>

	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list