jmsigner / amt

37 stars 13 forks source link

Simple shapefile conversion #22

Closed coverton-usgs closed 4 years ago

coverton-usgs commented 4 years ago

I apologize for what I am sure is a simple question but is there a convenient workflow for converting a track to shapefile with metadata intact? I can convert to other formats and eventually get to a writeOGR but I also run into a lot of errors that way:

example_sp<-as_ltraj(example, infolocs=example)
example_sp<-ltraj2spdf(example_sp)

gives Error in `[.data.frame`(tr, !is.na(tr$x), c("x", "y")) : undefined columns selected

but

example_move<-as_move(example)
example_sp<-move2ade(example_move)
example_sp@data<-example

works for making a spdf

jmsigner commented 4 years ago

Sorry, for the delay. I did not get notified for some reason. I am not sure I am sure what you want to do. Do you want to convert a track to a Shapefile with lines?

coverton-usgs commented 4 years ago

No worries at all I appreciate that you're willing to devote time to it at all, I know how things go!

Yes, a conversion from track to shapefile (either as points or lines, they both have their use)

jmsigner commented 4 years ago

I added to function to the development version of the package. You can source them for time being:

library(amt)
source("https://raw.githubusercontent.com/jmsigner/amt/master/R/coercion.R")

data(deer)

deer %>% as_sf_points() %>% sf::st_write("~/tmp/pts.shp")
deer %>% as_sf_lines() %>% sf::st_write("~/tmp/lines.shp")

Does this do what you are after?

coverton-usgs commented 4 years ago

The as_sf_points is working with an warning message:

Warning message:
In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  GDAL Message 6: Field t_ create as date field, though DateTime requested.

But the as_sf_lines is giving the error:

 Error in MtrxSet(x, dim, type = "MULTILINESTRING", needClosed = FALSE) : 
  is.list(x) is not TRUE 

Both rdgal and sf were older version on my machine but updating them had no effect on either the warning or error. Traceback below:

15.
stop(simpleError(msg, call = if (p <- sys.parent(1L)) sys.call(p))) 
14.
stopifnot(is.list(x)) 
13.
MtrxSet(x, dim, type = "MULTILINESTRING", needClosed = FALSE) 
12.
sf::st_multilinestring(l) at coercion.R#111
11.
sf::st_sfc(sf::st_multilinestring(l)) at coercion.R#111
10.
sf::st_sf(sf::st_sfc(sf::st_multilinestring(l))) at coercion.R#111
9.
as_sf_lines.track_xy(.) at coercion.R#88
8.
as_sf_lines(.) 
7.
function_list[[i]](value) 
6.
freduce(value, `_function_list`) 
5.
`_fseq`(`_lhs`) 
4.
eval(quote(`_fseq`(`_lhs`)), env, env) 
3.
eval(quote(`_fseq`(`_lhs`)), env, env) 
2.
withVisible(eval(quote(`_fseq`(`_lhs`)), env, env)) 
1.
z %>% as_sf_lines() %>% sf::st_write("D:/zlines.shp") 
jmsigner commented 4 years ago

Did you get the same error with the example data set? Could you email me or post the output of str(z)?

coverton-usgs commented 4 years ago

Yes, that would seem to be the problem. The example set works fine. My data has a character column for "id" instead of a numeric "burst" column.

Classes ‘track_xyt’, ‘track_xy’, ‘tbl_df’, ‘tbl’ and 'data.frame':  777109 obs. of  4 variables:
 $ x_: num  -122 -122 -122 -122 -122 ...
 $ y_: num  38.1 38.1 38.1 38.1 38.1 ...
 $ t_: POSIXct, format: "2015-01-20 17:06:47" "2015-01-20 17:08:22" "2015-01-20 17:10:09" "2015-01-20 17:11:56" ...
 $ id: chr  "X48511914641.WATE.01.1" "X48511914641.WATE.01.1" "X48511914641.WATE.01.1" "X48511914641.WATE.01.1" ...
 - attr(*, "crs_")=Formal class 'CRS' [package "sp"] with 1 slot
  .. ..@ projargs: chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"
coverton-usgs commented 4 years ago

However, converting the id to a numeric still gives me the same error:

z2<- z[,c("x_","y_","t_")]
z2$burst<- as.numeric(factor(z$id))
z2 %>% as_sf_lines() %>% sf::st_write("D:/z2lines.shp")

Neither changing projection nor limiting to a single individual resolves the issue either

jmsigner commented 4 years ago

This should be fixed now. Thanks for reporting back.

coverton-usgs commented 4 years ago

Thank you very much for the help, it does work now. The line file doesn't contain attributes, point file does, and I am having difficulty adding them using merge (its duplicating each record for every burst or id, but I think I am getting closer. Many thanks!

jmsigner commented 4 years ago

The problem is, that we would have to 'collapse' the attributes some how, which I think is difficult. Or which kind of attributes do you envision there?

coverton-usgs commented 4 years ago

This is really clunky (I'm a pretty poor coder), but this is what I was looking for:


linesegment<-function(in.file=x, id = "burst_", x_="x_", y_="y_", input_crs= get_crs(deer), outputlocation=getwd()) {
  x<-in.file
  x$RowID<-seq(from=1, to=nrow(x))
  DID<-as.data.frame(table(x[,id]))
  DID <- subset(DID, Freq > 0)
  DID$id<-as.character(DID$Var1)

  for(t in 1:nrow(DID)){
    Duck.temp<-x[x[,id]==DID[t,"id"],]
    if (nrow(Duck.temp)<= 1) next
    long <- as.list(Duck.temp[,x_])
    lat <- as.list(Duck.temp[,y_])
    coords <- matrix(c(long, lat), ncol = 2)
    linecoord<-data.frame(Longfrom=numeric(0), Latfrom=numeric(0),Longto=numeric(0), Latto=numeric(0),RowID=integer(0))
    for(j in 2:nrow(Duck.temp)){
      k=j-1
      linecoord[k, "Longfrom"] <- Duck.temp[k,"x_"]
      linecoord[k, "Latfrom"] <- Duck.temp[k,"y_"]
      linecoord[k, "Longto"] <- Duck.temp[j,"x_"]
      linecoord[k, "Latto"] <- Duck.temp[j,"y_"]
      linecoord[k,"RowID"]<- Duck.temp[j,"RowID"]
    }
    From<-data.frame(long=linecoord$Longfrom,lat=linecoord$Latfrom)
    To<-data.frame(long=linecoord$Longto, lat=linecoord$Latto)
    l<-vector("list", nrow(linecoord))
    for (i in seq_along(l)) {
      l[[i]] <- Lines(list(Line(rbind(From[i, ], To[i,]))), linecoord[i,"RowID"])
    }
    lsp<-SpatialLines(l, proj4string = input_crs)
    lspdf<-SpatialLinesDataFrame(lsp, data=Duck.temp[-1,], match.ID = FALSE)
    writeOGR(lspdf, dsn=outputlocation, layer=paste0(deparse(substitute(in.file)),"_",DID[t,3],"_ID_Line"), driver="ESRI Shapefile", overwrite_layer = TRUE)
    if(t==1){lspdf_comb<-lspdf}
    if(t>1){lspdf_comb<-bind(lspdf_comb,lspdf)}
        }
  writeOGR(lspdf_comb, dsn=outputlocation, layer=paste0(deparse(substitute(in.file)),"_combined_line_2"), driver="ESRI Shapefile", overwrite_layer = TRUE)
}

linesegment(in.file = deer,id = "burst_", x_="x_", y_="y_", input_crs = get_crs(deer), outputlocation = "D:/")
jmsigner commented 4 years ago

Cool, thanks a lot for sharing. I will have a look at it and try to include this into the package, if OK with you. But it may take some time.

coverton-usgs commented 4 years ago

Absolutely, Since I worked on this in my official capacity I should be getting approvals internally before providing it, but if it will take you a little while then I don't see any harm with planning on it. I did add an option to bypass the individual shapefile creation for each unique id and will endeavor to add some better error controls in the function as well. As soon as I get official word I'll send a push notification, otherwise feel free to adapt and incorporate as needed.

Thank you again very much for your time and effort on this package, I am starting to use it more and more.

Cory

coverton-usgs commented 4 years ago

Hello Johannes,

I just wanted to let you know that I received provisional approvals for release. I am hoping that they have made this repo public by now, If I make any changes to the code it will be reflected here:

https://code.usgs.gov/wildlife-telemetry-toolkit/telemetry-code-and-collaborations

Many thanks for your help, C