[STATSGRASS] Re: [GRASSLIST:10737] krigging result looks weird
Carlos "Guâno" Grohmann
carlos.grohmann at gmail.com
Tue Mar 7 19:17:34 EST 2006
Something like:
coords <- coordinates(srtm)
jcoords <- cbind(jitter(coords[,1]), jitter(coords[,2]))
> str(jcoords)
num [1:60888, 1:2] 744797 744896 744981 745063 745148 ...
> class(jcoords)
[1] "matrix"
> jcoords
[,1] [,2]
[1,] 744796.6 7277799
[2,] 744896.2 7277797
[3,] 744981.4 7277802
[4,] 745063.3 7277800
[5,] 745147.9 7277805
[6,] 745233.0 7277802
[7,] 745350.3 7277807
[8,] 745436.6 7277799
[9,] 745506.7 7277812
[10,] 745604.9 7277810
[11,] 745706.2 7277793
[12,] 745787.0 7277796
[13,] 745882.2 7277820
[14,] 745965.8 7277800
[15,] 746057.1 7277790
[16,] 746154.6 7277810
[17,] 746227.9 7277817
[18,] 746321.9 7277812
[19,] 746423.8 7277820
[20,] 746490.9 7277799
[21,] 746593.1 7277792
[22,] 746684.7 7277811
[23,] 746762.5 7277794
[24,] 746871.6 7277811
[25,] 746965.6 7277814
[26,] 747040.6 7277807
[27,] 747143.3 7277789
[28,] 747215.3 7277819
[29,] 747310.1 7277800
[30,] 747416.2 7277813
So they won't fall in the exac spot as before:
> OK_pred <- krige(srtm$cat ~ 1, loc=~jcoords[,1]+jcoords[,2], newdata=mask_SG, model=vrg.eye2, maxdist=270)
Error in "coordinates<-"(`*tmp*`, value = ~jcoords[, 1] + jcoords[, 2]) :
no direct or inherited method for function 'coordinates<-' for this call
I need a data frame:
> cat<-as.data.frame(srtm$cat)
> class(cat)
[1] "data.frame"
> srtm2<-SpatialPointsDataFrame(jcoords,cat)
> class(srtm2)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
> str(srtm2)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
..@ data :Formal class 'AttributeList' [package "sp"] with 1 slots
.. .. ..@ att:List of 1
.. .. .. ..$ srtm$cat: num [1:60888] 103 114 123 93 81 75 68 75 89 105 ...
..@ coords.nrs : num(0)
..@ coords : num [1:60888, 1:2] 744797 744896 744981 745063 745148 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : NULL
.. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
..@ bbox : num [1:2, 1:2] 744779 7261949 775681 7277821
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. ..@ projargs: chr NA
> OK_pred <- krige(srtm$cat ~ 1, loc=srtm2, newdata=mask_SG, model=vrg.eye2, maxdist=270)
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim, :
data items in gstat object and newdata have different
coordinate reference systems
I fotgot the proj4 when I created srtm2...
>srtm2<-SpatialPointsDataFrame(jcoords,cat,proj4string=CRS(G$proj4));
> str(srtm2)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
..@ data :Formal class 'AttributeList' [package "sp"] with 1 slots
.. .. ..@ att:List of 1
.. .. .. ..$ srtm$cat: num [1:60888] 103 114 123 93 81 75 68 75 89 105 ...
..@ coords.nrs : num(0)
..@ coords : num [1:60888, 1:2] 744798 744885 744968 745067 745171 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : NULL
.. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
..@ bbox : num [1:2, 1:2] 744779 7261949 775681 7277821
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. ..@ projargs: chr "+proj=utm +south +zone=22 +a=6378160
+rf=298.25 +no_defs +towgs84=-60.0,-2.0,-41.0"
Seems good, let's try it.
> OK_pred <- krige(srtm$cat ~ 1, loc=srtm2, newdata=mask_SG, model=vrg.eye2, maxdist=270)
[using ordinary kriging]
Hurray!!! This worked!
see the images:
www.igc.usp.br/pessoais/guano/tempo/krigging_good.jpg
www.igc.usp.br/pessoais/guano/tempo/krigging_good_zoom.jpg
Thanks Roger and Edzer!
Carlos
On 3/7/06, Carlos Guâno Grohmann <carlos.grohmann at gmail.com> wrote:
> And Yes, the artifacts are present in R.
>
> Carlos
>
> On 3/7/06, Carlos Guâno Grohmann <carlos.grohmann at gmail.com> wrote:
> > OK I did understant what you wan tme to do and why, but I didn't get how...
> > BTW I overlaid the vector points and the krigged map in GRASS and yes,
> > the artifacs position match the points position
> >
> > Carlos
> >
> > On 3/7/06, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > > On Tue, 7 Mar 2006, Edzer J. Pebesma wrote:
> > >
> > > > At the risk of being too obvious/ignorant, is the issue that kriging
> > > > reproduces
> > > > observed values where grid cell centres exactly coincide with data points?
> > > > A solution to that may be to shift the grid (or data) over a very small
> > > > distance.
> > >
> > > No, I asked Carlos to move this here because I suspected that there might
> > > be a direct explanation. Carlos, could you try jittering the input data
> > > points, would just moving them within a target 30m cell be enough, say
> > > within +/-15m? x+runif(*, -15, +15), y the same, and hope that both
> > > runif()s don't say 0?
> > >
> > > Thanks!
> > >
> > > Roger
> > >
> > > > --
> > > > Edzer
> > > >
> > >
> > > --
> > > Roger Bivand
> > > Economic Geography Section, Department of Economics, Norwegian School of
> > > Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> > > Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> > > e-mail: Roger.Bivand at nhh.no
> > >
> > >
> >
> >
> > --
> > +-----------------------------------------------------------+
> > Carlos Henrique Grohmann - Guano
> > Geologist M.Sc - Doctorate Student at IGc-USP - Brazil
> > Linux User #89721 - carlos dot grohmann at gmail dot com
> > +-----------------------------------------------------------+
> > _________________
> > "Good morning, doctors. I have taken the liberty of removing Windows
> > 95 from my hard drive."
> > --The winning entry in a "What were HAL's first words" contest judged
> > by 2001: A SPACE ODYSSEY creator Arthur C. Clarke
> >
>
>
> --
> +-----------------------------------------------------------+
> Carlos Henrique Grohmann - Guano
> Geologist M.Sc - Doctorate Student at IGc-USP - Brazil
> Linux User #89721 - carlos dot grohmann at gmail dot com
> +-----------------------------------------------------------+
> _________________
> "Good morning, doctors. I have taken the liberty of removing Windows
> 95 from my hard drive."
> --The winning entry in a "What were HAL's first words" contest judged
> by 2001: A SPACE ODYSSEY creator Arthur C. Clarke
>
--
+-----------------------------------------------------------+
Carlos Henrique Grohmann - Guano
Geologist M.Sc - Doctorate Student at IGc-USP - Brazil
Linux User #89721 - carlos dot grohmann at gmail dot com
+-----------------------------------------------------------+
_________________
"Good morning, doctors. I have taken the liberty of removing Windows
95 from my hard drive."
--The winning entry in a "What were HAL's first words" contest judged
by 2001: A SPACE ODYSSEY creator Arthur C. Clarke
More information about the grass-stats
mailing list