Skip to content

kde output from ks reversed

12 messages · andrewaduff, Robert J. Hijmans

Original
andrewaduff
# 
Im having trouble with the results from the ks package being reversed.  There
must be a simple error that I am overlooking?  Does anybody see my problem

## for testing, just the first 100 records, format x coord, y coord
curset <- curset[1:100,] 

#Calculate hpi, hpidiag
Hpi1 <- Hpi(x = curset)

## Do the kde
ud_hpi <- kde(x = curset, H = Hpi1)

im <- list(x=ud_hpi$eval.points[[1]], y=ud_hpi$eval.points[[2]],
z=ud_hpiv$estimate)
SPim <- image2Grid(im)

writeAsciiGrid(SPim,
"P:/GISProjects/Support_Projects/JimWatson/GOEA_ForJim_20110201/BBMM/TestBBMM/R/Final_KS/ud_hpi_v2.asc",
attr=ncol(SPim))


http://r-sig-geo.2731867.n2.nabble.com/file/n6037090/ReverseCapture.jpg
Robert J. Hijmans
# 
Perhaps I could have seen what you are overlooking if you had provided a
simple self-contained example. That should be very easy to make in this
case. It is not helpful to start an example script with 

curset <- curset[1:100,] 

This can be rather confusing as you never really know if the raster is
stored from top to bottom, and/or rotated. You could try if this works (I
believe it does):

library(raster)
r <- raster(ud_hpi)
plot(r)
r <- writeRaster(r, 'test.asc')
andrewaduff
# 
Here is a self contained sample with a different set of points, same issues
though.  When I use your suggested routines for export, I get a slightly
different issue but one I was seeing before my original post.  I avoided
this issue by using the code from my original post create the raster
manually,  Somehow the raster is getting warped in x and y???  When I use
your routines they fail, saying that the cells unequal in horizontal and
vertical dimensions.

*****************************

library(ks)
library(raster)
curset <- read.delim("sample.csv", header = TRUE, sep = ",")

#Calculate hpi
Hpi1 <- Hpi(x = curset, pre="scale")
Hpi1

## Do the kde
ud_hpi <- kde(x = curset, H=Hpi1)

library(raster)
r <- raster(ud_hpi)
plot(r)
writeRaster(r, 'test_r.asc')

*****************************
upon writeRaster I get an error - 

"Error in .startAsciiWriting(raster, filename, ...) : 
  x has unequal horizontal and vertical resolutions. Such data cannot be
stored in arc-ascii format"

http://r-sig-geo.2731867.n2.nabble.com/file/n6041112/sample.csv sample.csv 
http://r-sig-geo.2731867.n2.nabble.com/file/n6041112/sample.r sample.r
andrewaduff
# 
The more I think about it, I feel like something is going awry with the
coordinate system of the points I am using and the routines with respect to
scaling in hpi or kde.  I can see that the points in R are not where they
are supposed to be when I plot the UD  as such

plot(ud_hpi, drawpoints=TRUE, drawlabels=FALSE, col=3, lwd=3, cex=0.1)

If the plot points to the ArcGIS representation of the points it is wrong.

Those points are formatted in NAD83 Albers...

http://r-sig-geo.2731867.n2.nabble.com/file/n6041277/CoordinateCapture.jpg
Robert J. Hijmans
# 
though.
Thanks for trying, but this is not self-contained, as it refers to a file
that only you have.
That is a limitation of the ascii format, you can write to a different
format. 


Here is a self contained example that *looks* OK to me. 

library(ks) 
library(raster) 
set.seed=1
curset <- cbind(rnorm(100), rnorm(100))
Hpi1 <- Hpi(x = curset, pre="scale") 
ud_hpi <- kde(x = curset, H=Hpi1) 
r <- raster(ud_hpi) 
plot(r) 
points(curset)
writeRaster(r, 'test_r.tif')
andrewaduff
# 
"Thanks for trying, but this is not self-contained, as it refers to a file
that only you have. "

Im not sure what the problem is Robert, my previous post included both the
code file and my albers coordinates.  I will attach them again. Simulated
data doesn't really prove anything to me if you can't apply the results back
to the earth.

http://r-sig-geo.2731867.n2.nabble.com/file/n6041928/sample.csv sample.csv 

http://r-sig-geo.2731867.n2.nabble.com/file/n6041928/sample.r sample.r
Robert J. Hijmans
# 
code file
Got them now
? really?

This seems to work fine. Run your script and do 
plot(r)
points(curset)

You can create a proj.4 string from the info in your .prj file and do
projection(r) <- 'valid proj.4 string'
Or write r to a file and set the projection in ArcGIS.

r <- writeRaster(r, 'file.tif')

you cannot use ascii because the cells are not square.

Robert
andrewaduff
# 
Got it.  Thanks Robert.  What are good raster formats and R libraries for
dealing with non-square cells?  I guess im not used to working with them
much.

With respect to the UD values and subsequent GIS analysis - Lots of resource
utilization papers have been taking the hpi values from R and apply them in
Hawthorne Beyers 'Hawth's tool's' which writes to a square cell format. Is
this proper use of the hpi values from R?  How can I effectively get these
rasters back to square cells for comparitive home range analysis (VI, BA)
and resource selection with other 30 m data without altering the structure
and values of the UD?

Thanks for all your help.
Robert J. Hijmans
# 
cells?

You can use geotiff, as in my example ("file.tif").
I do not know that
I can see two ways to make the cells square. 

1) resampling the result
 
r2 = raster(r)
res(r2) = min(res(r))
r2 <- resample(r, r2, method='bilinear')


2) A better solution is to provide equally spaced evaluation points, using a
pre-defined square-cell grid (here re-using r2):

eval <- xyFromCell(r2, 1:ncell(r2))
ud_hpi <- kde(x = curset, H=Hpi1, eval.points=eval) 
r3 <- setValues(r2, ud_hpi$estimate)


Robert
andrewaduff
# 
Thanks for all your help Robert - I appreciate it as Im sure you a busy man!
Ill check with the author of that tool as it seems his conversions most
logically would have been using your resampling suggestion in ArcObjects. 
Since I am looking at a number of methods for home range analysis including
brownian bridge, fixed kernels, and others.  I want to make sure that I am
comparing the resultant surfaces correctly.
This seems less than optimal, but as long as the values were rescaled to a
true pdf it seems that it could work fine.

2) A better solution is to provide equally spaced evaluation points, using a
pre-defined square-cell grid 

I do have concern that while all three rasters (r, r2, and r3) show the same
spatial pattern, the values of the sum of the surface are different.  I
guess the key here is rescaling by the sum of the cell values so that they
sum to 1 (a true pdf) and then I can create a volume surface.

 ## Get the volume
udvol_hpi <- getvolumeUD(r)

It looks like I will have to work on getting these rasters into class
"khrud" or "kbbhrud".  Ill do some research on that and if I figure it out I
will post.  All these different software packages (within R and outside R)
can make home range comparison tricky.  Adehabitat goes a long ways to bring
everything together.