r-spatial / lwgeom

bindings to the liblwgeom library
https://r-spatial.github.io/lwgeom/
58 stars 23 forks source link

Add st_startpoint #22

Closed Robinlovelace closed 6 years ago

Robinlovelace commented 6 years ago

I think this has the high-level R code in place - just need to figure out how to make it work on the C/Rcpp side of things. Using this as a learning opportunity as much as anything else and new to cpp so apologies if some of the changes here don't make sense (I know the 4.1 in lwgeom.cpp makes no sense - I'm experimenting to see different error messages as I'll comment below).

Help on how to fix this appreciated - heads-up @virgesmith.

Robinlovelace commented 6 years ago

Heads-up @edzer and @dbaston this seems to be working now as illustrated in the reprex below. Next step I think: convert the matrix back to sf object but that should be relatively straightforward.

devtools::install_github("robinlovelace/lwgeom", ref = "startpoint")
#> Using GitHub PAT from envvar GITHUB_PAT
#> Skipping install of 'lwgeom' from a github remote, the SHA1 (48a022a6) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(sf)                                                         
#> Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2
m = matrix(c(0, 1, 2, 0, 1, 4), ncol = 2)                           
l = st_sfc(st_linestring(m))                                        
lwgeom::st_startpoint(l)                                            
#>      [,1] [,2]
#> [1,]    0    0
l2 = st_sfc(st_linestring(m), st_linestring(m[3:1, ]))              
lwgeom::st_startpoint(l2)                                           
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    2    4
Session info ``` r devtools::session_info() #> Session info ------------------------------------------------------------- #> setting value #> version R version 3.5.0 (2018-04-23) #> system x86_64, linux-gnu #> ui X11 #> language en_GB:en #> collate en_GB.UTF-8 #> tz Europe/London #> date 2018-06-04 #> Packages ----------------------------------------------------------------- #> package * version date source #> backports 1.1.2 2017-12-13 CRAN (R 3.5.0) #> base * 3.5.0 2018-04-23 local #> class 7.3-14 2015-08-30 cran (@7.3-14) #> classInt 0.2-3 2018-04-16 cran (@0.2-3) #> compiler 3.5.0 2018-04-23 local #> curl 3.2 2018-03-28 cran (@3.2) #> datasets * 3.5.0 2018-04-23 local #> DBI 1.0.0 2018-05-02 cran (@1.0.0) #> devtools 1.13.5 2018-02-18 CRAN (R 3.5.0) #> digest 0.6.15 2018-01-28 CRAN (R 3.5.0) #> e1071 1.6-8 2017-02-02 cran (@1.6-8) #> evaluate 0.10 2016-10-11 CRAN (R 3.3.2) #> git2r 0.21.0 2018-01-04 CRAN (R 3.5.0) #> graphics * 3.5.0 2018-04-23 local #> grDevices * 3.5.0 2018-04-23 local #> grid 3.5.0 2018-04-23 local #> htmltools 0.3.6 2017-04-28 CRAN (R 3.5.0) #> httr 1.3.1 2017-08-20 cran (@1.3.1) #> knitr 1.20 2018-02-20 CRAN (R 3.5.0) #> lwgeom 0.1-5 2018-06-04 Github (robinlovelace/lwgeom@48a022a) #> magrittr 1.5 2014-11-22 CRAN (R 3.3.2) #> memoise 1.1.0 2017-04-21 CRAN (R 3.5.0) #> methods * 3.5.0 2018-04-23 local #> R6 2.2.2 2017-06-17 cran (@2.2.2) #> Rcpp 0.12.17 2018-05-18 CRAN (R 3.5.0) #> rmarkdown 1.9 2018-03-01 CRAN (R 3.5.0) #> rprojroot 1.3-2 2018-01-03 CRAN (R 3.5.0) #> sf * 0.6-3 2018-05-17 cran (@0.6-3) #> spData 0.2.8.9 2018-06-03 Github (nowosad/spData@7b53933) #> spDataLarge 0.2.6.5 2018-06-02 Github (nowosad/spDataLarge@bc058ad) #> stats * 3.5.0 2018-04-23 local #> stringi 1.2.2 2018-05-02 CRAN (R 3.5.0) #> stringr 1.3.1 2018-05-10 CRAN (R 3.5.0) #> tools 3.5.0 2018-04-23 local #> udunits2 0.13 2016-11-17 cran (@0.13) #> units 0.5-1 2018-01-08 cran (@0.5-1) #> utils * 3.5.0 2018-04-23 local #> withr 2.1.2 2018-03-15 CRAN (R 3.5.0) #> yaml 2.1.19 2018-05-01 CRAN (R 3.5.0) ```
Robinlovelace commented 6 years ago

One question @edzer would you suggest returning the result as an sfc object? Not sure how st_startpoint() would deal with multi-dimensional points which seems to be the purpose of the POINT4D class. Feedback appreciated - will be useful in the stplanr pkg is the motivation.

edzer commented 6 years ago

That seems more useful than a matrix, yes. The logic is to replace the geometry if the object is of class sf. I don't know how to deal with Z and/or M coordinates - sth to sort out!

Robinlovelace commented 6 years ago

OK great. Do you have any example data representing Z/M input data? Should add that to the @examples section also I imagine.

dbaston commented 6 years ago

Next step I think: convert the matrix back to sf object but that should be relatively straightforward.

These LWPOINT constructors may be useful: https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/liblwgeom.h.in#L1374

edzer commented 6 years ago
sf::read_sf(system.file("shape/storms_xyz_feature.shp", package = "sf"))

There is also an xyzm there.

cmcaine commented 6 years ago

Robin asked me to take a look at this. I'm glad this is on its way in - I had to write my own R function for this problem last time I used sf (converting the line into points and using head wasn't performant).

I don't know how to deal with Z and/or M coordinates

Probably everyone worked this out already or this obvious solution misses the point, but in case not:

This function takes an sfc object and they know how many dimensions they have, so this function can just ask the input what dimensions it has and ignore the response of lwgeom/startpoint for dimensions that will not be meaningful. Similarly, the dimensions of the input tell us what LWPOINT constructor to use.

As a fun aside, lwgeom gives NO_Z_VALUE or NO_MVALUE results for z and m if the given geometry doesn't have enough dimensions. NO*_VALUE = NO_VALUE = 0.0 and a decade old comment suggests that it will change to NaN at some point!

Robinlovelace commented 6 years ago

Thanks for the input! Interesting to read a 10 year old comment from Paul Ramsey. In terms of next steps I suggest we check it out on Tuesday, unless there's an easy fix in there.