SpatialLinesDataFrame to owin in spatstat (by Dr. Anderson)

Have you ever wondered how to convert a SpatialLinesDataFrame to an object of type owin in ‘spatstat’? In most cases, it would not make sense to covert a line file to a window object, but there may be some situations (as described below) when you may need to do it. As far as I know, however, there is no existing function in the ‘spatstat’ library that can accomplish this task, but with a little backdoor ingenuity, it can be done.

I figured this out (the hard way) when I was trying to import a boundary file to use as analytical window for a point pattern analysis. After using ‘rgdal’ to import a shapefile (via readOGR()), I was able to read in the  shapefile, but I was confused when it would not convert to an object of type owin in ‘spatat’ (via the standard owin() or as.owin()). When I checked the class of the object I imported with readOGR, it turned out to be a SpatialLinesDataFrame object rather than a SpatialPolygonsDataFrame object. I figured this must be the problem.

LLP <- readOGR(dsn=getwd(), layer=”LLP_boundary”)
OGR data source with driver: ESRI Shapefile
Source: “C:\Users\coreanderson\Documents”, layer: “LLP_boundary”
with 1 features
It has 2 fields

> as.owin(LLP)
Error in as.owin.default(Z, …, fatal = fatal) :
Can’t interpret W as a window

> class(LLP)
[1] “SpatialLinesDataFrame”
attr(,”package”)
[1] “sp”

When I edited the GPS file in ArcGIS, I must have saved it as a polyline shapefile rather than a polygon shapefile. This is an easy fix in ArcGIS, but I was working from home (with no access to the license server), so ArcGIS was not an option. I decided to crack into the SpatialLinesDataFrame object using str() to see if I could pull out the coordinates to build the owin object from scratch.

> str(LLP)
Formal class ‘SpatialLinesDataFrame’ [package “sp”] with 4 slots
..@ data :’data.frame’: 1 obs. of 2 variables:
.. ..$ LEFT_FID : int -1
.. ..$ RIGHT_FID: int 0
..@ 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:395, 1:2] 284003 284002 284002 284002 284002 …
.. .. .. ..@ ID : chr “0”
..@ bbox : num [1:2, 1:2] 283778 3401613 284228 3402002
.. ..- 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=17 +datum=WGS84 +units=m +no_defs +ellps=WG

While I fully appreciate ‘sp’ as a foundational R package for conducting spatial analysis, accessing a slot in an object from ‘sp’ is not intuitive. Above you can see that there are two slots with coordinate data, one is a sublist in the slot @lines and the other is a sublist in the slot @bbox; it is (presumably) the former that we want. The question is: how do you pull it out? The only advise I can offer is that it is useful to remember that S4 objects have slots that act similarly to lists and to keep tinkering until you figure it out. After some trial and error, I finally cracked the code:

> LLP@lines[[1]]@Lines[[1]]@coords

This produced a matrix with two columns representing the X and Y coordinates (in UTMs). In inspecting the matrix, I noticed that the first and last coordinates were the same (which is the standard format for a polygon); however, ‘spatstat’ specifies that the first and last coordinates should be different. I decided to isolate the X and Y columns and then remove last line in each.

> XX <- LLP@lines[[1]]@Lines[[1]]@coords[,1] 
> XX <- XX[-395]
> YY <- LLP@lines[[1]]@Lines[[1]]@coords[,2] 
> YY <- YY[-395]

Now that I had the coordinates isolated, I thought I could build the owin object, but I was not quite there (of course).

> Z <- owin(poly=list(x=XX, y=YY))
Error in owin(poly = list(x = XX, y = YY)) :
Area of polygon is negative – maybe traversed in wrong direction?

Luckily, ‘spatstat’ catches the exemption here and throws an informative error message. Now, if I can only reverse those vectors…

> Z <- owin(poly=list(x=rev(XX), y=rev(YY)))  #Bingo!; plot(z)

I constantly remind my students that writing code is a “two-steps forward one-step back process”, but with a little logic and steadfast grit, most things are possible. The next step is to write a little function to automate this process; stay tuned.

Leave a Comment

Filed under Uncategorized

Leave a Reply

Your email address will not be published. Required fields are marked *