centroids and masking UDs

[Thomas Collins FAL]

Hi all, 

I’d like to add my solution to a new thread concerning both the masking of UDs to impassable features (E.g. rivers described here and here ) and taking measurements of the centre of UDs. This is because firstly it was painful to complete but also because I believe there are questions linking these things together. 

My study species is found on a peninsula. I am interested in the centre of the home range as well movement away during key times. I want to measure distance between a centroid against the locations of key events. For this, I wanted to mask to the shape of the peninsula as there is overlap into water. I had thought once masked, I would get a new measurement of the mode and mean. However, and of course!, the mean and the mode described here are the same after masking. I am now considering my original central measurement which was the centre of 50% polygon.

The painful part of masking was due to partly to flipping axes as well as a confusion with ‘cols’ and ‘rows’ in my ‘UD’:

# so… to mask UD to spdf…

# change the projection of spdf to ‘UD’:

mweya_sp <- sp::spTransform(mwey_sp, CRSobj = ctmm::projection(UD))

 

# make new raster same col/row as UD:

rr <- raster::raster(ncol = 195, nrow = 149)

 

# get  extent of ‘UD’

xmin = raster::extent(UD)[1,1]

xmax = raster::extent(UD)[2,1]

ymin = raster::extent(UD)[1,2]

ymax = raster::extent(UD)[2,2]

 

# change extent of ‘rr’ to extent of ‘UD’:

raster::extent(rr) <- c(xmin, xmax, ymin, ymax)   

 

# rasterize ‘mweya-sp’ to ‘rr’ :

rr <- raster::rasterize(mweya_sp, rr, field = 1)

 

# flip it (BECAUSE OTHERWISE IT WILL NOT MASK THE CORRECT AREAS AKA UPSIDE DOWN)!

rr <- raster::flip(rr, direction = 'y')

 

# resample to ‘UD’ to align to ‘UD’:

rr_re <- raster::resample(rr, ctmm::raster(UD))

 

# Change values so that we have ‘1’ for ‘shape’ and ‘NA’ for ‘not-shape’:

rr_re[rr_re][rr_re[rr_re] > 0] <- 1

 

# for loop to mask to ‘shape’:

  for(i in 1:ncol(UD$CDF)) {

    for(j in 1:nrow(UD$CDF)) {

       if(is.na(raster::as.matrix(rr_re)[i,j])) {

                    UD$CDF[j, i] <- 1                          # i and j are reversed here

      }

    }

  }  # this works, hooray !!

 

But.... to check the cols and rows:

ncol(UD$CDF)

  y

  195

 

nrow(UD$CDF)

  X

149

 

Compare this with the output of the raster of the ‘UD’:

ctmm::raster(UD)

 

class      : RasterLayer

dimensions : 195, 149, 29055  (nrow, ncol, ncell)      *** here cols and rows are reversed

resolution : 21.76689, 36.70524  (x, y)

extent     : -1214.177, 2029.09, -2656.904, 4500.617  (xmin, xmax, ymin, ymax)

crs        : +proj=tpeqd +lat_1=-0.182323355226995 +lon_1=29.8912322690628 +lat_2=-0.184633847485593 +lon_2=29.8937692160309 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

source     : memory

names      : layer

values     : 0.001060171, 1  (min, max)

 

I am unsure why I had to flip my axes and reverse I and j for the for loop but my feeling is that it is because of projection problems (or a problem with rookie coding!). Does anyone an explanation for this? 

 

Lastly, and to add to the conversation about centroid measurements  - if I am masking the UD, which should I use to determine the centre of the akde output; mean, mode or the centroid of the 50% polygon (*or other)? 

 Best wishes and thanks for this amazingly helpful resource!

Tom 

 

 

 

[Christen Fleming]

Hi Tom,

Another way to do this would be to ask if (x=UD$r[i],y=UD$r[j]) is in the polygon (after transforming it to the same projection as the UD).

In general, I think of the mode as being the center of the home range. This would not change much from masking (slightly if the masking was performed at the kernel level).

Best,

Chris

[Thomas Collins FAL]

Hi Chris,

thanks for your thoughts there. After finding that my raster loop wasn't lining up correctly wit han spdf i tried this loop for an spdf object of the peninsula:

 for(i in 1:nrow(UD$CDF)) {       

          for(j in 1:ncol(UD$CDF)) {     

            point <- data.frame(x = UD$r$x[i], y = UD$r$y[j])   

            coordinates(point) <- c("x", "y")

            proj4string(point) <- proj4string(mweya_sp_PR)

            T_F <- rgeos::gContains(mweya_sp_PR, point)    

            

            if((T_F) == F) {

              UD$CDF[i, j] <- 1                         

            }

          }

        }

its ever so slow! any ideas how to improve the speed of this loop?

best,

Tom

[Christen Fleming]

Hi Tom,

I suspect the slowest part of the calculation is rgeos::gContains(). You might try some competing functions and seeing which are faster, like sp::point.in.polygon().

If this was code that you were going to run a million times, there are also some mild speed ups that you can get from vectorizing instead of looping, by constructing the data.frame first with something like outer() and then feeding the entire thing to the evaluation function, and then taking the output of that to mask the CDF.

Best,

Chris

[Thomas Collins FAL]

thanks for that. ill give it a go.

best,

Tom