2

I hope you can help me with this problem i can't find how to overcome. Sorry if I made some mistakes while writing this post, my english is a bit rusty right now.

Here is the question. I have .shp data that I want to analyze in R. The .shp can be either lines that represent lines of traps we set to catch octopuses or points located directly over those lines, representing where we had catured one.

Clarification image

The question i'm trying to answer is: Are octopuses statistically grouped or not?

After a bit of investigation it seems to me that i need to use R and its linearK function to answer that question, using the libraries Maptools, SpatStat and Sp.

Here is the code i'm using in RStudio:

Loading the libraries

library(spatstat)
library(maptools)
library(sp)

Creating a linnet object with the track

t1<- as.linnet(readShapeSpatial("./20170518/t1.shp"))

I get the following warning but it seems to work

Warning messages:
1: use rgdal::readOGR or sf::st_read 
2: use rgdal::readOGR or sf::st_read

Plotting it to be sure everything is ok

plot(t1)

Plot of t1 object

Creating a ppp object with the points

p1<- as.ppp(readShapeSpatial("./20170518/p1.shp"))

I get the same warning here, but the real problems start when I try to plot it:

> plot(p1)
Error in if (!is.vector(xrange) || length(xrange) != 2 || xrange[2L] <  : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: Interpretation of arguments maxsize and markscale has changed (in spatstat version 1.37-0 and later). Size of a circle is now measured by its diameter. 
2: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) :
  All mark values are NA; plotting locations only.
3: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) :
  All mark values are NA; plotting locations only.
4: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) :
  All mark values are NA; plotting locations only.
5: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) :
  All mark values are NA; plotting locations only.
6: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) :
  All mark values are NA; plotting locations only.
7: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) :
  All mark values are NA; plotting locations only.

Now what is left is to join the objects in a lpp object and to analyze it with the linearK function

> pt1 <- lpp(p1,t1)
> linearK(pt1)
Function value object (class ‘fv’)
for the function r -> K[L](r)
......................................
    Math.label     Description        
r   r              distance argument r
est {hat(K)[L]}(r) estimated K[L](r)  
......................................
Default plot formula:  .~r
where “.” stands for ‘est’
Recommended range of argument r: [0, 815.64]
Available range of argument r: [0, 815.64]

This is my situation right now. What i dont know is why the plot function is not working with my ppp object and how to understant the return of the linearK function. Help(linearK) didn't provide any clue. Since i have a lot of tracks, each with its set of points, my desired outcome would be some kind of summary like x tracks analized, a grouped, b dispersed and c unkown.

Thank you for your time, i'll greatly appreciate if you can help me solve this problem.

Edit: Here is a link to a zip file containing al the shp files of one day, both tracks and points, and a txt file with my code. https://drive.google.com/open?id=0B0uvwT-2l4A5ODJpOTdCekIxWUU

Bhargav Rao
  • 50,140
  • 28
  • 121
  • 140
Alberto Mnemon
  • 117
  • 4
  • 14

2 Answers2

1

First two pieces of general advice: (1) each time you create a complicated object, print it at the terminal, to see if it is what you expected. (2) When you get an error, immediately type traceback() and copy the output. This will reveal exactly where the error is detected.

A ppp object must include a specification of the study region (window). In your code, the object p1 is created by converting data of class SpatialPointsDataFrame, which do not include a specification of the study region, converted via the function as.ppp.SpatialPointsDataFrame, into an object of class ppp in which the window is guessed by taking the bounding box of the coordinates. Unfortunately, in your example, there is only one data point in p1, so the default bounding box is a rectangle of width 0 and height 0. [This would have been revealed by printing p1.] Such objects can usually be handled by spatstat, but this particular object triggers a bug in the function plot.solist which expects windows to have non-zero size. I will fix the bug, but...

In your case, I suggest you do 

Window(p1) <- Window(t1)

immediately after creating p1. This will ensure that p1 has the window that you probably intended.

If all else fails, read the spatstat vignette on shapefiles...

Adrian Baddeley
  • 1,956
  • 1
  • 8
  • 7
  • Thank you so much Adrian for your answer. I have also noticed that the ppp object is where all my troubles begin, so I decided to try creating it in a more manual way, and it seem to work. – Alberto Mnemon Aug 24 '17 at 10:04
0

I have managed to find a solution. As Adrian Baddeley noticed there was a problem with the owin object. That problem seems to be bypassed (not really solved) if I create the ppp object in a manual way instead of converting my set of points.

I have also changed the readShapeFile function for the rgdal::readOGR, since the first once was deprecated, and that was the reason of the warnings I was getting.

This is the R script i'm using right now, commented to clarify:

#first install spatstat, maptools y sp

#load them

library(spatstat)
library(maptools)
library(sp)


#create an array of folders, will add more when everything works fine
folders=c("20170518")
for(f in folders){

  #read all shp from that folder, both points and tracks
  pointfiles <- list.files(paste("./",f,"/points", sep=""), pattern="*.shp$")
  trackfiles <- list.files(paste("./",f,"/tracks", sep=""), pattern="*.shp$")

  #for each point and track couple
  for(i in 1:length(pointfiles)){


    #create a linnet object with the track

    t<- as.linnet(rgdal::readOGR(paste("./",f,"/tracks/",trackfiles[i], sep="")))
    #plot(t)



    #create a ppp object for each set of points

    pre_p<-rgdal::readOGR(paste("./",f,"/points/",pointfiles[i], sep=""))
    #plot(p)
    #obtain the coordinates the current set of points
    c<-coordinates(pre_p)
    #create vector of x coords
    xc=c()
    #create vector of y coords
    yc=c()
    #not a very good way to fill my vectors but it works for my study area
    for(v in c){
      if(v>4000000){yc<-c(yc,v)}
      else {if(v<4000000 && v>700000){xc<-c(xc,v)}}
    }
    print(xc)
    print(yc)
    #create a ppp object using the vectors of x and y coords, and a window object
    #extracted from my set of points
    p=ppp(xc,yc,Window(as.ppp(pre_p)))

    #join them into an lpp object

    pt <- lpp(p,t)
    #plot(pt)

    #analize it with the linearK function, nsim=9 for testing purposes
    #envelope.lpp is the method for analyzing linear point patterns
    assign(paste("results",f,i,sep="_"),envelope.lpp(pt, nsim=9, fun=linearK))

  }#end for each points & track set
}#end for each day of study

So as you can see this script is testing for CSR each couple of points and track for each day, working fine right now. Unfortunately I have not managed to create a report or reportlike with the results yet (or even to fully understand them), I'll keep working on that. Of course I can use any advice you have, since this is my first try with R and many newie mistakes will happen.

The script and the shp files with the updated folder structure can be found here(113 KB size)

Alberto Mnemon
  • 117
  • 4
  • 14