University of California, Los AngelesDepartment of StatisticsStatistics C173/C273 Instructor: Nicolas ChristouSpatial data analysis with GRASS and R - ExampeWe will use the Wolfcamp data from Texas. First we need to create a new location that will holdthese data. You can name it wf, and also create a mapset a1 . Read the data first in R to find thenorth-south and east-west extend of the region. To create the location in GRASS we use projectionvalues and choose x, y coordinates. The minimum value of x is -145.2 and its maximum 112.8. Theminimum value of y is 9.4 and its maximum 184.8. Therefore define the location as follows:North: 186South: 8West: -144East: 114and give resolution 1 to both N to S and E to W.You are done!Enter GRASS through the location wf and mapset a1. Once you are at the GRASS prompt accessthe ascii file wolfcamp.txt that you have received through e-mail and you have saved on yourcomputer’s Desktop. To import this file you use the v.in.ascii command. Give an output namewolfcamp in. See the snapshots of the v.in.ascii module below:1Define the columns of the data set:Once you have imported the file, display the points (first open a monitor window, and then typed.vect wolfcamp in). You can customize the size and the type of the symbols as shown below.For more information type d.vect help. Here we used:GRASS 6.3.cvs (wf):~ > d.vect wolfcamp_in size=5 icon=basic/circleto get the following plot:2From GRASS to R:Now we want to take these data into R, analyze them using gstat, and then get the krigingpredictions again back to GRASS. Simply at the command line of GRASS we type R. Once you arein R load the libraries gstat and spgrass6 as follows:> library(gstat)> library(spgrass6)You can now read the data from GRASS into R by typing:a <- readVECT6("wolfcamp_in")If you type a at the R prompt you will see your data. The first 5 rows are shown here:> a[1:5,]coordinates cat x y level1 (42.8, 127.6, 0) 1 42.8 127.6 14642 (-27.4, 90.8, 0) 2 -27.4 90.8 25533 (-1.2, 84.9, 0) 3 -1.2 84.9 21584 (-18.6, 76.5, 0) 4 -18.6 76.5 24555 (96.5, 64.6, 0) 5 96.5 64.6 1756We observe that the data are read as “SpatialPointsDataFrame” (the number in parenthesis areread as coordinates).Now, we can compute variograms, fit variogram models, and do kriging predictions. For this exam-ple we use universal kriging (see handout on universal kriging). Here are some useful commands:g <- gstat(id="level", formula = level~x+y, data = a)q <- variogram(g)plot(q)v.fit <- fit.variogram(q, vgm(30000,"Sph",60,10000))plot(q, model=v.fit)> v.fitmodel psill range1 Nug 9959.434 0.000002 Sph 33852.314 66.08075Here is the fitted model variogram to the sample variogram:3We create now the grid for the kriging predictions:> x.range <- as.integer(range(a$x))> x.range[1] -145 112> y.range <- as.integer(range(a$y))> y.range[1] 9 184grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1),y=seq(from=y.range[1], to=y.range[2], by=1))class(grd)coordinates(grd) <- ~x+ygridded(grd) <- TRUEclass(grd)plot(grd, cex=0.5)points(a, pch=1, col="green", cex=0.9)The grid with the data points:4Now we update the gstat object and perform universal kriging as follows:> g <- gstat(g, id="level", model=v.fit) #Update object g.> pr <- predict.gstat(g, model=v.fit, newdata=grd)[using universal kriging]> names(pr)[1] "level.pred" "level.var"> pr$level.pred> pr$level.varFrom R back to GRASS:Now we can save the predicted values and their variances as follows:> writeRAST6(pr, "level_pred", zcol="level.pred")Projection of input dataset and current location appear to match.Proceeding with import...100%> writeRAST6(pr, "level_var", zcol="level.var")Projection of input dataset and current location appear to match.Proceeding with import...100%Type q() to leave R and go back to GRASS.In GRASS open a new monitor to display the raster maps.d.mon x0d.rast level_predd.rast level_varYou can also create a vector map that displays the contours with the following command:GRASS 6.3.cvs (wf):~ > r.contour level_pred out=level_pred_cont step=50Reading data:100%Range of data: min = 785.597107 max = 3647.292236Range of levels: min = 800.000000 max = 3600.000000Displacing data:100%Total levels: 57WARNING: 97 crossings foundsBuilding topology ...66 primitives registered5Building areas: 100%0 areas built0 isles builtAttaching islands:Attaching centroids: 100%Topology was built.Number of nodes : 66Number of primitives: 66Number of points : 0Number of lines : 66Number of boundaries: 0Number of centroids : 0Number of areas : 0Number of isles : 0r.contour complete.Here is the raster map with the contours and the data points:GRASS 6.3.cvs (wf):~ > d.rast level_pred100%GRASS 6.3.cvs (wf):~ > d.vect level_pred_contGRASS 6.3.cvs (wf):~ > d.vect wolfcamp_in size=5 icon=basic/circle6A very interesting and useful tool for a 3D visualization in GRASS is the nviz tool. At the GRASSprompt you typeGRASS 6.3.cvs (wf):~ > nviz level_predThis command will display the following:If you want to display the data points click on “Visualize” and access the data points through“Vector Points”. Also, using the options “perspective”, “z-exag”, and “height” we get the followingthree-dimensional view plot from the northeast corner.Using the available options under “Appearance” we can change the colors, add text to the plot,add fringes, etc. You can also save this display by clicking on “File” and then on “Save image as”and you can save it as a TIFF