[R-sig-Geo] Bearing angle of UTM projected SpatialLines
Eduardo Diez
eduardodiez at gmx.com
Wed May 25 21:16:34 CEST 2016
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
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
>>
>
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list