[R-sig-Geo] ODP: ODP: spsample and NA
Edzer J. Pebesma
e.pebesma at geo.uu.nl
Mon Jun 25 11:45:36 CEST 2007
Not only to avoid NAs, but also to get approximately right sample size.
Stratified sampling does not take place based on strata defined by
attributes (map layers), but rather by dividing the area into square
blocks (strata) and randomly sample (n=1) from each block. That is why
the resulting sample doesn't have exactly the prescribed/requested size
(although one could iterate). See also Ripley 1981 Spatial Statistics
for explanation.
Best wishes,
--
Edzer
Roger Bivand wrote:
> There are two forms of sp objects for grids, SpatialPixelsDataFrame and SpatialGridDataFrame, and you need the former to avoid NAs. I am not sure how to proceed without access to a copz of zour data to see whz spsample is including the NA cells.
>
> Roger
>
> --- Roger Bivand, NHH, Helleveien 30, N-5045 Bergen, Roger.Bivand at nhh.no
>
> ________________________________
>
> Od: Agustin Lobo w imieniu Agustin Lobo
> Wys³ano: Cz 21.06.2007 09:41
> Do: Roger Bivand
> DW: R-sig-Geo at stat.math.ethz.ch
> Temat: Re: ODP: [R-sig-Geo] spsample and NA
>
>
>
> The test with meuse.grid works:
>
> #copied from file:SpatialGridDataFrame-class.html
> coordinates(meuse.grid) = c("x", "y") # promote toSpatialPointsDataFrame
> gridded(meuse.grid) <- TRUE # promote to SpatialPixelsDataFrame
> x = as(meuse.grid, "SpatialGridDataFrame") # creates the full grid
> #Now me:
> > unique(x at data$soil)
> [1] NA 1 2 3
>
> > table(x at data$soil)
> 1 2 3
> 1665 1084 354
>
> test.stsamp <- spsample(x, n = 10, "stratified")
> o <- overlay(x, test.stsamp)
> > o
> coordinates part.a part.b dist soil ffreq
> 5556 (179156, 330917) 0 1 0.1030270 1 1
> 6670 (180039, 330336) 0 1 0.0703333 1 2
> 2240 (180677, 332606) 1 0 0.0975136 1 1
>
>
> But don't know how to make that my raster behave in the same way. What
> does it mean "cast to SpatialPixels first"? I thought that as
> X comes from an imported raster it was an sp object already.
>
> Also, there are several components in the x coming from meuse.grid:
> str(x)
> Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
> ..@ data :'data.frame': 8112 obs. of 5 variables:
> .. ..$ part.a: num [1:8112] NA NA NA NA NA NA NA NA NA NA ...
> .. ..$ part.b: num [1:8112] NA NA NA NA NA NA NA NA NA NA ...
> .. ..$ dist : num [1:8112] NA NA NA NA NA NA NA NA NA NA ...
> .. ..$ soil : num [1:8112] NA NA NA NA NA NA NA NA NA NA ...
> .. ..$ ffreq : Factor w/ 3 levels "1","2","3": NA NA NA NA NA NA NA
> NA NA NA ...
>
> how does spsample know which layer has to be used for the stratification?
>
> Thanks!
>
> Agus
>
> Roger Bivand escribió:
>
>> I think that if you cast to SpatialPixels first, the problem will go away, because the all-NA cases are dropped with their coordinates. Could you try the same example with meuse.grid and see how it goes?
>>
>> Roger
>>
>> ---
>> Roger Bivand, NHH, Helleveien 30, N-5045 Bergen, Roger.Bivand at nhh.no
>>
>> ________________________________
>>
>> Od: r-sig-geo-bounces at stat.math.ethz.ch w imieniu Agustin Lobo
>> Wys³ano: ¦r 20.06.2007 10:11
>> Do: R-sig-Geo at stat.math.ethz.ch
>> Temat: [R-sig-Geo] spsample and NA
>>
>>
>>
>> It is often the case, at least with raster (grid) objects,
>> that you are not interested on a part of the raster.
>> For example because it's outside a non.rectangular region of study
>> and/or because some types of terrain (some classes) are outside
>> the scope of the study. In any case, those parts can be
>> represented by NA.
>>
>> Unfortunately, while some functions like table() would not
>> include the NA in the statistics, spsample includes
>> the NA as just another class, which implies that a fraction
>> of the resulting coordinates will lay on the NA
>> part of the raster. This makes it difficult, to get N valid (i.e., not
>> on NA parts) points using spsample. Is there any way around this?
>> For example, if I want 10 points on non-NA cells of
>> raster X, I would like
>>
>> spsample(X, 10, "stratified")
>>
>> to provide 10 positions on non-NA cells.
>>
>> At the moment, I get:
>>
>> x <- GDAL.open("C:/ALOBO/dipu2006/STLL_VEG2006/ESTRAT3_OBAC.tif")
>> ESTRAT3OBAC <- asSGDF_GROD(x, output.dim=c(100,100))
>> GDAL.close(x)
>>
>> ESTRAT3OBAC at data[ESTRAT3OBAC at data==55537]<- NA
>>
>> X <- ESTRAT3OBAC
>> X$cut <- as.ordered(cut(X$band1, c(0,10,100),include.lowest=TRUE))
>> table(X$cut)
>>
>> [0,10] (10,100]
>> 83 1961
>>
>> summary(X$cut)
>> [0,10] (10,100] NA's
>> 83 1961 7956
>> table(X$cut)
>>
>> [0,10] (10,100]
>> 83 1961
>>
>> Xspsample <- spsample(X,10,"stratified")
>> Xo <- overlay(X,Xspsample)
>> summary(Xo)
>> Object of class SpatialPointsDataFrame
>> Coordinates:
>> min max
>> x1 410995.6 414954.4
>> x2 4608002.1 4610566.0
>> Is projected: TRUE
>> proj4string : [ +proj=utm +zone=31 +ellps=intl +units=m +no_defs]
>> Number of points: 9
>> Data attributes:
>> band1 cut
>> Min. : 67 [0,10] :0
>> 1st Qu.: 97 (10,100]:3
>> Median : 159 NA's :6
>> Mean :24745
>> 3rd Qu.:55537
>> Max. :55537
>>
>> That is, I get 6 positions on NA values of X
>>
>> Thanks!
>>
>> Agus
>>
>> --
>> Dr. Agustin Lobo
>> Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
>> LLuis Sole Sabaris s/n
>> 08028 Barcelona
>> Spain
>> Tel. 34 934095410
>> Fax. 34 934110012
>> email: Agustin.Lobo at ija.csic.es
>> http://www.ija.csic.es/gt/obster
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>
> --
> Dr. Agustin Lobo
> Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
> LLuis Sole Sabaris s/n
> 08028 Barcelona
> Spain
> Tel. 34 934095410
> Fax. 34 934110012
> email: Agustin.Lobo at ija.csic.es
> http://www.ija.csic.es/gt/obster
>
>
>
> [[alternative HTML version deleted]]
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list