TransWikia.com

0........0.kriging error messages, grid file problem ? (automap, gstat)

Geographic Information Systems Asked by Maria Eugenia D'Amato on January 15, 2021

I am trying to conduct kriging of variable called percent in South Africa.

   head(LL12)

     Location Longitude  Latitude   percent
    1       L1  16.48476 -28.60737 1.0000000
    2       L3  16.90045 -29.24154 1.0000000
    3       L4  16.99268 -28.44906 0.6363636
    4       L5  17.06871 -29.67825 1.0000000
    5       L6  17.09935 -29.00169 0.6000000
    6       L7  17.25630 -28.82418 1.0000000

I got the grid as follows:

       SAgrd <- makegrid(SApol_border, n = 1000)
       colnames(SAgrd) <- c('lon','lat')

    SAgrd_pts <- SpatialPoints(coords = SAgrd, 
                       proj4string=CRS(proj4string(SApol_border)))


    SAgrd_pts_in <- SAgrd_pts[SApol_border, ]
    plot(SAgrd_pts_in)

enter image description here

# change to shorter name and transform 
        SA.grid = SAgrd_pts_in

# transform grid to s.pix.d.f.
        SA.grid=as(SA.grid,"SpatialPixelsDataFrame") 

# Ordinary kriging, no new_data object (automap)
# obtain a spdf

        coordinates(LL12) <- ~ Longitude + Latitude
        LL12spdf <- LL12  # indicate the data is a spdf

        krigLL12 = autoKrige(percent~1, LL12spdf)

thus far, good. The above renders a graphic.
enter image description here

However I want to obtain a better graphic where the points are within the frame boundaries/shape of the country. Then I fit a variogram with automap:

    LL12vgm <- autofitVariogram(percent~1, LL12spdf) 

Then use the fitted variogram and the grid in ordinary kriging (https://rdrr.io/cran/automap/man/autoKrige.html):

    L12_am_krg = autoKrige(percent~1, LL12spdf, SA.grid)

however I get the following error message:

Error in autoKrige(percent ~ 1, LL12spdf, SA.grid) : 
  Either input_data or new_data is in LongLat, please reproject.
   input_data:  NA 
   new_data:    +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
In addition: Warning messages:
1: In proj4string(new_data) :
  CRS object has comment, which is lost in output
2: In proj4string(new_data) :
  CRS object has comment, which is lost in output
3: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
 but +towgs84= values preserved

Then I tried using gstat and the fitted variogram :

      L12.kriged <- krige(percent ~ 1, LL12spdf, SA.grid, model=LL12vgm$var_model)

and I get the following error message
1 NA
1 "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"

Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
  var1 : data item in gstat object and newdata have different coordinate reference systems
In addition: Warning messages:
1: In proj4string(newdata) :
  CRS object has comment, which is lost in output
2: In proj4string(newdata) :
  CRS object has comment, which is lost in output

I would prefer to use the automap for its simplicity.
It looks like I need to reformat the grid file but I cannot figure out how.
I am a novice in these topics.

One Answer

I managed to get both krige (gstat) and autoKrige (automap) to work erasing the ref sysyem form the grid file and getting a SpatialPixels file as follows:

    SAg <- makegrid(SApol_border, n = 1000)
    colnames(SAg) <- c('Longitude','Latitude')

    SAg_pts <- SpatialPoints(coords = SAg, 
                        proj4string=CRS(proj4string(SApol_border)))

    SAg_pts_in <- SAg_pts[SApol_border, ]
    class(SAg_pts_in) # Spatial Points
    crs(SAg_pts_in) # check coord ref system
    crs(SAg_pts_in) <- NA
    class(SAg_pts_in)  # [1] "SpatialPoints"
                               attr(,"package")
                         [1] "sp"

    SAgrd <- as(SAg_pts_in, "SpatialPixels") 

    L12_am_krg = autoKrige(percent~1, LL12spdf, SAgrd)
    plot(L12_am_krg)

enter image description here

    L12.kriged <- krige(percent ~ 1, LL12spdf, SAgrd, model=LL12vgm$var_model)
    plot(L12.kriged)

enter image description here

Answered by Maria Eugenia D'Amato on January 15, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP