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

Paul Hiemstra p.hiemstra at geo.uu.nl
Thu Dec 6 21:28:29 CET 2007


Hi Dale,

Going from latlong to utm destroys the grid yeah. I used the inverse 
distance weighted interpolation of the gstat package. Along the lines of:

# 1. Create target grid, for example by using a bbox polygon of your
# area of interest using spsample. Getting a bbox-polygon from a 
SpatialObject
# little home brewn function :)

bbox2polygon = function(obj) {
# This function converts a bounding box of a spatial object to a SpatialPolygons-object
# obj :  Spatial object
  bb = bbox(obj)
  pl = Polygon(data.frame(rbind(c(bb["x","min"], bb["y","max"]),
                                c(bb["x","max"], bb["y","max"]),
                                c(bb["x","max"], bb["y","min"]),
                                c(bb["x","min"], bb["y","min"]),
                                c(bb["x","min"], bb["y","max"])
                            )))
  pls = Polygons(list(pl), ID = "bbox")
  spat_poly = SpatialPolygons(list(pls), proj4string = CRS(proj4string(obj)))
  return(spat_poly)
}


data(meuse)
coordinates(meuse) =~x+y
bb = bbox2polygon(meuse)
spplot(meuse, "zinc",sp.layout= list("sp.polygons",bb))
grd = spsample(bb, cellsize = c(50,50), type = "regular") # Cell size in m
spplot(meuse, "zinc",sp.layout= list("sp.points",grd))

# 2. Interpolate to this grid, mind to keep either the number of max
# allowd points small (nmax), or limit the search radius (maxdist). 
Experiment a bit
require(gstat)
new_grd = krige(zinc~1, meuse, grd, maxdist = 50) # In your case you 
ofcourse
    # have much more points and more regularly spaced
gridded(new_grd) = TRUE
spplot(new_grd)

Good luck!

cheers,
Paul

Dale Steele schreef:
>> gridded(Depth)=TRUE
>>     
> suggested tolerance minimum: 0.323243088554591Error in
> points2grid(points, tolerance) :
>   dimension 1 : coordinate intervals are not constant
>
> I suspect this implies that my 'Depth' data are not located on a
> regular grid.   I originally downloaded from the internet, and
> converted from long/lat to UTM using PBSMapping.  Do I need to perform
> some sort of interpolation?
>
> Best Regards ,
>     --Dale
>
>
> On Dec 6, 2007 2:11 PM, Paul Hiemstra <p.hiemstra at geo.uu.nl> wrote:
>   
>> Hi Dale,
>>
>> I think you need to change the SpatialPointsDataFrame Depth to a
>> SpatialGrid. Use:
>>
>> gridded(Depth) = TRUE
>>
>>
>> cheers,
>> Paul
>>
>> Dale Steele schreef:
>>     
>>> 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
>>>>
>>>>
>>>>
>>>>         
>> --
>>
>> 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
>>
>>
>>     


-- 
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