[R-sig-Geo] How to find ocean depth at specific locations?

Dale Steele dale.w.steele at gmail.com
Thu Dec 6 20:07:45 CET 2007


Thank you very much.  I think I've made progress, but am now stuck....

I now have two data frames (sample appended below).  I would like to
find the indices of 'Depth' on locations in 'Scallop' so I can add
'depth' as a covariate at each sample location in 'Scallop'.

Two issues arise when I try code adapted from your example.

1) The code below generates an error message, when I try to find the
indices of 'Depth' at locations in 'Scallop'.

> overlay(Depth, Scallop)
Error in function (classes, fdef, mtable)  :
  unable to find an inherited method for function "overlay", for
signature "SpatialPointsDataFrame", "SpatialPointsDataFrame"

2) The code below first plots the points from Scallop, then draws
Depth (color image) which covers(obscures) the points drawn
previously.

spplot(Depth, "depth", sp.layout= list("sp.points", Scallop))
  -
2) The code below generates an error message, when I try to find the
indices of 'Depth' at locations in 'Scallop'.

> overlay(Depth, Scallop)
Error in function (classes, fdef, mtable)  :
  unable to find an inherited method for function "overlay", for
signature "SpatialPointsDataFrame", "SpatialPointsDataFrame"

>Depth[1:5,]
          coordinates depth
7  (593.212, 4539.22)    10
22 (614.239, 4539.51)     4
23 (615.643, 4539.53)     6
24 (617.048, 4539.56)     8
25 (618.444, 4539.58)     8

> Scallop[1:5,]
         coordinates tcatchp1
1 (792.142, 4494.53)        1
2 (795.331, 4485.39)        1
3 (778.164, 4490.29)        1
4 (767.395, 4475.07)        2
5 (773.325, 4467.87)        1

Thanks.  --Dale

On Dec 2, 2007 4:26 PM, Paul Hiemstra <p.hiemstra at geo.uu.nl> wrote:
> Hi Dale,
>
> If you have all your data in Spatial objects (sp-package) and they are
> both in the same projection you can use the overlay function. Read
> ?overlay for more info. Small example:
>
> # Script to extract the distance to the river meuse (holland) from a
> grid (meuse.grid) for the points in meuse.
> library(sp)
> data(meuse)
> coordinates(meuse) =~x+y
> data(meuse.grid)
> gridded(meuse.grid) = ~x+y
> spplot(meuse.grid, "dist", sp.layout = list("sp.points", meuse))
> # Get indices of meuse.grid that contain the values of the grid at the
> overlay(meuse.grid, meuse)
> # Get the values from the grid and attach to meuse
> meuse$new_dist = meuse.grid$dist[overlay(meuse.grid, meuse)]
>
> hope this helps!
>
> cheers,
> Paul
>
> Dale Steele schreef:
>
> > I would like to use ocean depth as a covariate in a spatial model
> > where I have a measured catch at a specific location, but depth at
> > that site was not reported in the original dataset.  The dataset
> > 'scallop' in package SemiPar is a representative example:
> >
> >
> >> scallop[20:25,]
> >>
> >    latitude longitude tot.catch
> > 20 39.46667 -72.66667       722
> > 21 39.43333 -72.73333        41
> > 22 39.40000 -72.80000        39
> > 23 39.31667 -72.73333         0
> > 24 39.33333 -72.85000        94
> > 25 39.26667 -72.81667         7
> >
> >
> > I've been able to download a very large grid of depths in this area
> > (33514 points) from: <http://topex.ucsd.edu/cgi-bin/get_data.cgi>.
> > Given this data, I'm stuck at how to estimate the depth at each
> > location in the scallops dataset.
> >
> > depths <- read.table("depths.txt", header=F,col.names=c("x","y","z"))
> > depths$x <- depths$x - 360
> >    ## subset the data keeping only negative depths
> > depths <- subset(depths, depths$z < 0)
> >    ## Convert depths to a positive number
> > depths$z <- (-1)*depths$z
> > length(depths$z)
> >
> > Thanks.  --Dale
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
> --
> Drs. Paul Hiemstra
> Department of Physical Geography
> Faculty of Geosciences
> University of Utrecht
> Heidelberglaan 2
> P.O. Box 80.115
> 3508 TC Utrecht
> Phone:     +31302535773
> Fax:    +31302531145
> http://intamap.geo.uu.nl/~paul
>
>




More information about the R-sig-Geo mailing list