[R-sig-Geo] ODP: ODP: spsample and NA
Agustin Lobo
Agustin.Lobo at ija.csic.es
Mon Jun 25 13:02:15 CEST 2007
I understand, but how do I convert between one type and the other?
I've tried:
x <- GDAL.open("C:/ALOBO/dipu2006/STLL_VEG2006/ESTRAT3_OBAC.tif")
ESTRAT3OBAC <- asSGDF_GROD(x, output.dim=c(118,189))
GDAL.close(x)
ESTRAT3OBAC at data[ESTRAT3OBAC at data==55537]<- NA
X <- ESTRAT3OBAC[seq(1,118,by=10),seq(1,189,by=10)] #small example
> summary(X)
Object of class SpatialGridDataFrame
Coordinates:
min max
x 410649.7 415399.7
y 4607910.5 4610910.5
Is projected: TRUE
proj4string : [ +proj=utm +zone=31 +ellps=intl +units=m +no_defs]
Number of points: 2
Grid attributes:
cellcentre.offset cellsize cells.dim
x 410774.7 250 19
y 4608035.5 250 12
Data attributes:
band1
Min. : 4.0
1st Qu.: 67.0
Median :156.0
Mean :125.8
3rd Qu.:163.0
Max. :237.0
NA's :119.0
> Y <- SpatialPixels(X)
> summary(Y)
Object of class SpatialPixels
Coordinates:
min max
x 410649.7 415399.7
y 4607910.5 4610910.5
Is projected: TRUE
proj4string : [ +proj=utm +zone=31 +ellps=intl +units=m +no_defs]
Number of points: 2
So X and Y are not the same at all.
Thanks for your help.
Agus
Edzer J. Pebesma escribió:
> 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
>>
>
>
--
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
More information about the R-sig-Geo
mailing list