GRASS_R_meetup[1]

GRASS_R_meetup[1] - University of California, Los Angeles...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: University of California, Los Angeles Department of Statistics Introduction to GRASS/R spatial data analysis What is Geographic Information Systems (GIS)? - "A geographic information system (GIS) integrates hardware, software, and data for capturing, managing, analyzing, and displaying all forms of geographically referenced information." (ESRI, http://www.gis.com). - "GIS are tools that allow for the processing of spatial data into information and used to make decisions about some portion of the earth." (Fundamentals of Geographic Information Systems, DeMers. M., Third Edition, 2005) - "GIS is often described as integration of data, hardware, and software designed for management, processing, analysis, and visualization of georeferenced data." (OPEN SOURCE GIS, A GRASS GIS APPROACH, Neteler, M., Mitasova, H., Third Edition (2008). Why GIS? Needed for map analysis operations in areas such as - Forestry - Fire Departments - Police Departments - EMS - Traffic - Water resources - Real estate - Military - Business - Mining, geology - Hydrology - Wildlife management - Epidemiology - Agriculture - Election results - Etc. Useful tool for statisticians that work with spatial data. Geographical Data: Four main components. - Input - Administration - Analysis - Presentation 1 With GIS we can study - Land - Elevation - Forests - Population densities - Energy consumption - Mineral resources, etc. A GIS map has layers. Each layer represents geographic objects that are similar. For example, cities can be one layer, rivers another one, roads, lakes, etc. With a click of a button, we can include all these layers on our GIS map, or simply have as many layers as we want. GIS can help us understood our geographic data easily. 3D visualization and animation (see snapshot below from GRASS): 2 Spatial data - Geostatistics Why spatial statistics? Noel Cressie (in his textbook "Statistics for Spatial Data") writes "why, how, when" are not enough to describe our data. We need to add "where" (where did they occur?). First demonstration of spatial data appear in the form of maps (Halley (1686)). Spatial models appeared much later (Student (1907), R. A. Fisher (1935), Fairfield (1938), Whittle (1954), Matheron (1960s)), etc. Today, spatial statistics models appear in areas such as mining, geology, hydrology, ecology, environmental science, medicine, image processing, crop science, epidemiology, forestry, atmospheric science, etc. Need to develop models that deal with data collected at different spatial locations. The basic components are the spatial locations {s1 , s2 , , sn } and the data observed at these locations denoted as {Z(s1 ), Z(s2 ), , Z(sn )}. The distance between the observations is important in analyzing spatial data. With distance we mostly mean "Euclidean distance". However there are other forms of distances (e.g. road miles, travel time, etc.). The latter is modeled through multidimensional scaling. Here we will consider mostly (if not always) Euclidean distances. Geostatistics was developed for geology and mining, but in general it means statistics applied to problems in the earth (geo=earth, in Greek, ) sciences. A very popular method for analyzing spatial data is kriging, which was developed by Georges Matheron (1960s and 1970s). He coined the term kriging in honor of a South African mining engineer (Danie Krige). Terminology in kriging: Variogram, nugget, range, sill, correlogram, stationarity, variance-covariance matrix, spherical variogram, exponential variogram, ordinary kriging, universal kriging, block kriging, isotropic and anisotropic variogram, spatial prediction. Other methods of spatial prediction and interpolation. Raster map of the predicted values qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqq qqqq q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqq q qqqqq q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqq qqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqq q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q qq qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqq qq q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqq q q qqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q qq qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqq qqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqq qq qqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqq qq q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqq qqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqq q qq q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q qqq q q q qq q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqq qq q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqq qqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q qq q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q qqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q qq q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqq q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqq qqq qq q q qq qqqqqqqqqqqqqqqqqq q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqq qqq q qqq q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q qqq q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq qqq q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q qqqqq qqqqqqqqqqq q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q qq q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq 42 q 5 0.06 0.0 40 6 q q q q q q q q qq q 0.07 q q q q q qq South to North 0.05 q qq q q qq q q q qq q q q q q q 38 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.07 qq q q q q q qq 36 q q 0q0 . q q q q 5 q qq 0qq .0 q 8 q qq q q qq q q q 0.075 q q q q q 34 qq qq q q q 0. q q06 q qq q qq q qq q q .07 q q q 0 32 -124 -122 -120 West to East -118 -116 -114 32 q q 34 q q qq qq 0.07 q qqqqq qq q q q q qqqq q q q q q q q q q q q q 0.08 q 36 q q 0.08 q q q y q q q q qq 38 q 40 42 0.05 5 -124 -122 -120 x -118 -116 -114 3 General GIS principles It is important to understand the basic terminology and functionality of GIS (certain principles are common to all systems). Some widely used GIS softwares (commercial or open source) can be found at http://en.wikipedia.org/wiki/List_of_GIS_software http://www.freegis.org/ Georeferenced data include spatial component and attribute component. Spatial component: Two basic models. Field representation (raster data model): Things such as temperature, elevation, air pollution, vegetation, rain fall, wind speed do not have shapes. Instead they have measurable values for any particular location on the surface of earth. The most common surface is a raster. It is a matrix of equally sized square cells (pixels), where each cell represents a unit of surface area (for example, 1 square meter). A raster map that represents elevation may cover an area of 100 square kilometers. If there were 100 cells in this raster, each cell would represent one square kilometer of equal width and height (that is, 1km 1km). Each cell contains a measured or estimated value for that location. Here are two examples. Elevation: Snow level: Another type of raster maps are the so called thematic maps. Each cell gets a categorical value (for example, if we analyze vegetation data, 1, 2, 3, 4, 5 may represent different categories of vegetation). Here is one example. 4 The resolution of the raster map controls the level of spatial detail. The smaller the cell size the smoother the raster will be, however it takes longer to process it and requires larger storage space. On the other hand, if the cell size is too large, information may be lost. Geometrical objects representation (vector data model:): Geographic features are represented by lines (e.g. roads, rivers), points (e.g. cities, data points), polygons (e.g. county or state boundaries). This representation is given by the coordinates. The vector data model is based on arc-node representation, where an arc is a series of points given by (x, y) coordinates. The two endpoints of an arc are called nodes, while points along a line are called vertices. Two consecutive points define an arc segment. Here is one example. Can you guess what this vector map represents? Attributes: Attributes are descriptive data that provide information associated with raster or vector data. 5 Introduction to Geographical Resources Analysis Support System (GRASS) General information: GRASS Website: http://grass.fbk.eu/. GRASS logo: GRASS Team: http://grass.ibiblio.org/community/team.php. GRASS Wiki site: http://grass.osgeo.org/wiki/. GRASS manual: http://grass.ibiblio.org/gdp/html_grass64/ GRASS was originally developed for the U.S. Army Construction Engineering Research Laboratories (1982-1995) as a tool for land managing. GRASS is a free software used for data management, image processing, graphics production, spatial modelling, and visualization of many types of data. It is currently used in Academia, Industry, and Government (NASA, USGS, U.S. Census Bureau, etc.). The GRASS database: All data and files are saved under the home directory or under a shared network, e.g. ~/Users/nicolas/grass Under this directory, the GRASS GIS data are organized by projects stored in subdirectories called LOCATIONs. Each location is defined by its coordinates and map projection. Once a new location is created, GRASS internally produces other directories for the handling of the data (not important at this point). Further, LOCATIONs are subdivided into map subdirectories called MAPSETs. Note: When creating a new location GRASS automatically creates a special mapset called PERMANENT. Each LOCATION can have several MAPSETs. We can represent the above as follows: GRASS Database (Home directory) LOCATION /B MAPSET /PERMANENT/ b1 /A /C /PERMANENT /a1 /a2 /b2 /PERMANENT /c1 /c2 Under a mapset we can save our files that contain maps (either vector or raster with their attributes and other features). 6 GRASS modules (function classes): The modules in GRASS are organized by name (based on their function class). The first letter refers to the function class, then it follows by a dot and then by another word that describes the specific task described by the module. For example d.rast will display a raster map. Here are the most important modules in GRASS: First letter d. db. g. i. m. ps. r. r3. v. Function class display database general imagery misc postscript raster 3D raster vector Type of command graphical display database management general file operations imagery processing miscellaneous commands creation of a Postscript map raster data processing 3D raster data processing vector data processing The general syntax of a GRASS command is similar to UNIX commands. To get help in GRASS type the module and then "help". For example: d.rast help To open a monitor for map display we type: d.mon x0 Note: GRASS can open up to 7 monitors. If more than one monitors are opened, we can select a particular monitor as follows: d.mon select=x0 "A major pleasure in working with spatial data is their visualization. Maps are amongst the most compelling graphics, because the space they map is the space we think we live in, and maps can show things we cannot see otherwise". Roger Bivand, Edzer Pebesma, Virgilio Gomez-Rubio, Applied Spatial Data Analysis with R, Use R!, Springer 2008. 7 Examples of raster and vector data: Elevation map: A vector map: 8 Before we begin using the GRASS commnands we will first learn how to create a new LOCATION. There are three ways to create a new location (by clicking on Georeferenced file, or EPSG codes, or Projection values as shown on the screenshot below. We will begin using the last one (Projection values). 9 To create a new LOCATION using the last option, click on "Projection values". This will take you to the next window: Suppose we want to create the LOCATION "west1". Simply enter the name west1 next to LOCATION (as shown above) and hit return. Give a name to mapset (say, aaaaa). Then press ESC and then ENTER and follow the steps below as described by the screenshots... 10 11 Note: Here we choose the option "A" because we do not have any coordinates available at the moment. We choose a general non-georeferenced coordinate system x, y. 12 Note: Set the WEST EDGE and SOUTH EDGE values to zero, and enter for NORTH EDGE the number of rows and for EAST EDGE the number of columns. The next page explains what the number of columns and rows are. For the GRID RESOLUTION you can choose 1, because the units are pixels. See next page... 13 Note: To find the number of rows and columns: If you are using a MAC computer, click on Apple key+i. If you are using a PC computer, do right click and select properties. The example below shows the map of South America and the corresponding number of rows and columns. The number of rows and columns that we enter (see previous snapshots) should cover at least the size of the map that we want to import, It can be defined larger than needed. This will not use extra memory since the memory used depends only on the actual size of the file imported. 14 The map we want to import in GRASS for this tutorial is the following: 15 16 So far, we have created a new LOCATION (west1) and a new MAPSET (aaaaa). Also, GRASS has created automatically the PERMANENT mapset. See next snapshot of the GRASS startup screen. 17 We can now create new MAPSETs by first selecting the LOCATION we want (here, west1) and then entering the MAPSET's name in the box below "Create new mapset in selected location". Say, we want to create a new mapset named aa1. Here is the output: 18 To enter GRASS, we select the mapset "aaaaa" and click on "Enter GRASS". This will take you to the following: Note: The upper left window as shown in this snapshot, is the terminal "GRASS Command Window". This is the window where you type your commands (more will follow soon...). The upper right window it is called the "GIS Manager" and the window on the lower right is the "Map Display" window. The window on the lower left contains messages after we run a GRASS command. The last three windows are linked together. The "GRASS Command Window" can also display maps, etc., but we will have to open a monitor first. So, there are two ways to work with GRASS data: Either through the GRASS Command Window and a monitor, or with the other three windows as shown in the snapshot above. 19 Import and display raster data Suppose we want to import and display the map "snow.jpg", which located under ~/grass/snow.jpg On the terminal GRASS Window Command type the following commands (in the snapshots below, we can see the output of these commands): GRASS 6.3.cvs (w1):~ > r.in.gdal in=~/grass/snow.jpg out=snow Projection of input dataset and current location appear to match. Proceeding with import... 100% 100% 100% GRASS 6.3.cvs (w1):~ > g.list rast ---------------------------------------------<raster> files available in mapset <a3>: snow.blue snow.green snow.red ---------------------------------------------GRASS 6.3.cvs (w1):~ > d.mon x0 using default visual which is TrueColor ncolors: 16777216 Graphics driver [x0] started GRASS 6.3.cvs (w1):~ > d.rast snow.green 100% GRASS 6.3.cvs (w1):~ > r.composite GRASS 6.3.cvs (w1):~ > d.rast snow3 100% 20 Here is what each one of the previous GRASS commands does: > r.in.gdal in=~/grass/snow.jpg out=snow Imports the map snow.jpg and gives it the name "snow". > g.list rast Lists all the maps created by GRASS. Here we have snow.blue, snow.green, snow.red. > d.mon x0 Opens a new window in order to display the map. > d.rast snow.green Displays the map snow.green (see screenshot below). > r.composite Opens a new window in order to combine the three maps (see screenshot below). > d.rast snow3 Displays the map with the three colors combined. Map snow.green is displayed: Note: GRASS uses the RGB color model. In this model the three colors (red, green, blue) are used in varying intensities to produce other colors. The intensity of each of the three colors (red, blue, green) is measured on a scale from 0-255 (0 represents no color, 255 represents maximum intensity). For example, the black color can be obtained when R = 0, B = 0, G = 0, etc. 21 Combine snow.blue, snow.green, snow.red into a map named snow3: 22 The output: 23 The new map: > d.rast snow3 24 Another way to display maps is through the GIS Manager and the Map Display window. This is shown below in the 4 snapshots. First click on the raster icon below "Config". Then click on "Base map:" below and select the raster map snow3. 25 And looks like this: And you can display the map at the Map Display window (click on the upper left icon): 26 The g.region command: Managing map regions and resolutions To change the region, resolution, and boundaries of the map the g.region command is used with the following flags and options: -d, -p, -l, -e, -c, -o, -dp, res, n, s, w, e, save. -p -d -l -e -c save res= nsres= ewres= n= s= e= w= Print region Gives Default region Gives Print lat/long Gives extent of the region Gives center coordinates Saves the current region Changes resolution Changes resolution n-s Changes resolution e-w Changes the north extent Changes the south extent Changes the east extent Changes the west extent Here are some examples: Here is how we can get the current region numbers: GRASS 6.3.cvs (west1):~ > g.region -pec rast=snow3 projection: 0 (x,y) zone: 0 north: 1878 south: 0 west: 0 east: 900 nsres: 1 ewres: 1 rows: 1878 cols: 900 cells: 1690200 north-south extent: 1878.000000 east-west extent: 900.000000 center easting: 450.000000 center northing: 939.000000 GRASS 6.3.cvs (west1):~ > 27 Suppose we want to select a particular region: GRASS 6.3.cvs (west1):~ > g.region n=1000 e=500 rast=snow3 projection: 0 (x,y) zone: 0 north: 1000 south: 0 west: 0 east: 500 nsres: 1 ewres: 1 rows: 1000 cols: 500 cells: 500000 north-south extent: 1000.000000 east-west extent: 500.000000 center easting: 250.000000 center northing: 500.000000 GRASS 6.3.cvs (west1):~ > -pec We can go back to the default region as follows: g.region rast=snow3 or GRASS 6.3.cvs (west1):~ > g.region -dp projection: 0 (x,y) zone: 0 north: 1878 south: 0 west: 0 east: 900 nsres: 1 ewres: 1 rows: 1878 cols: 900 cells: 1690200 GRASS 6.3.cvs (west1):~ > 28 We can change the resolution of the region: GRASS 6.3.cvs (west1):~ > g.region nsres=2 ewres=2 rast=snow3 -pec projection: 0 (x,y) zone: 0 north: 1878 south: 0 west: 0 east: 900 nsres: 2 ewres: 2 rows: 939 cols: 450 cells: 422550 north-south extent: 1878.000000 east-west extent: 900.000000 center easting: 450.000000 center northing: 939.000000 GRASS 6.3.cvs (west1):~ > There are many options to set the region. Here is another one: The north and east extent are given in terms of the south and west. The result will be a region on the south-west corner of the current region. GRASS 6.3.cvs (west1):~ > g.region n=s+500 e=w+500 rast=snow3 -p projection: 0 (x,y) zone: 0 north: 500 south: 0 west: 500 east: 900 nsres: 1 ewres: 1 rows: 500 cols: 400 cells: 200000 north-south extent: 500.000000 east-west extent: 400.000000 center easting: 700.000000 center northing: 250.000000 GRASS 6.3.cvs (west1):~ > The changed region can be saved using g.region -s, but this must be done from the PERMANENT MAPSET. In the next pages we will explore how GRASS and the package gstat of R can be used to analyze spatial data. The method of spatial prediction called kriging will be used. The installation of the R/GRASS interface can be done by a single command: install.packages("spgrass6", "gstat", type="source", dependencies=TRUE) 29 Spatial data analysis with GRASS and R - Exampe We will use the data from the Walker Lake of Nevada (variable v - see figure below). The variable V is measured in parts per million (ppm) at the Walker Lake area in the state of Nevada (from An Introduction to Applied Geostatistics by Issaks, E., and Srivastava, R. M. (1989)). You can access the data at: a1 <- read.table("http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/ walker_lake_v.txt", header=TRUE) First we need to create a new location that will hold these data. You can name it w13, and also create a mapset a1 . Read the data first in R to find the north-south and east-west extend of the tegion. To create the location in GRASS we use projection values and choose x, y coordinates. The minimum value of x is 8 and its maximum 251. The minimum value of y is 8 and its maximum 291. Therefore define the location as follows: North: 252 South: 7 West: 7 East: 292 and give resolution 1 to both N to S and E to W. You are done! Enter GRASS through the location w13 and mapset a1. Once you are at the GRASS prompt the data can be accessed and imported using the v.in.ascii command. Give an output name v in. See the snapshots of the v.in.ascii module below: 30 Define the columns of the data set: Once you have imported the file, display the points (first open a monitor window, and then type d.vect v 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 (w13):~ > d.vect v_in size=5 icon=basic/circle to get the following plot: 31 From GRASS to R: Now we want to take these data into R, analyze them using gstat, and then get the kriging predictions again back to GRASS. Simply at the command line of GRASS we type R. Once you are in 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: walker_v <- readVECT6("v_in") If you type walker v at the R prompt you will see your data. The first 5 rows are shown here: > walker_v[1:5,] coordinates cat x y v 1 (11, 8, 0) 1 11 8 0.0 2 (8, 30, 0) 2 8 30 0.0 3 (9, 48, 0) 3 9 48 224.4 4 (8, 69, 0) 4 8 69 434.4 5 (9, 90, 0) 5 9 90 412.1 We observe that the data are read as "SpatialPointsDataFrame" (the number in parenthesis are read as coordinates). Now, we can compute the variogram, fit variogram models, and do kriging predictions. Here are some useful commands: g <- gstat(id="v", formula=v~1, data=walker_v) q <- variogram(g) plot(q) v.fit <- fit.variogram(q, vgm(60000,"Sph",60,30000)) plot(q, model=v.fit) Here is the fitted model variogram to the sample variogram: 32 We create now the grid for the kriging predictions: x.range x.range [1] 8 y.range y.range [1] 8 <- range(walker_v$x) 251 <- range(walker_v$y) 291 grd <- 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+y gridded(grd) <- TRUE class(grd) plot(grd, cex=0.5) points(walker_v, pch=1, col="green", cex=0.7) The grid with the data points: Now we update the gstat object and perform ordinary kriging as follows: g <- gstat(g, id="v", model=v.fit ) #Update object g. v_ok <- predict.gstat(g, model=v.fit, newdata=grd) #Or using the krige function: v_ok <- krige(id="v", v~1, walker_v, grd, model = v.fit) names(v_ok) v_ok$v.pred v_ok$v.var 33 From R back to GRASS: Now we can save the predicted values and their variances as follows: > writeRAST6(v_ok, "v_pred", zcol="v.pred") Projection of input dataset and current location appear to match. Proceeding with import... 100% > writeRAST6(v_ok, "v_var", zcol="v.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 maps. d.mon x0 d.rast v_pred d.rast v_var You can also create a vector map that displays the contours with the following command: GRASS 6.3.cvs (w13):~ > r.contour v_pred out=v_pred_cont step=50 Reading data: 100% Range of data: min = -78.295837 max = 1528.099976 Range of levels: min = 0.000000 max = 1500.000000 Displacing data: 100% Total levels: 31 WARNING: 85 crossings founds Building topology ... 652 primitives registered Building areas: 100% 0 areas built 0 isles built Attaching islands: Attaching centroids: 100% Topology was built. Number of nodes : 680 Number of primitives: 652 Number of points : 0 Number of lines : 652 Number of boundaries: 0 Number of centroids : 0 Number of areas : 0 Number of isles : 0 r.contour complete. Here is the raster map with the contours: GRASS 6.3.cvs (w13):~ > d.vect v_pred_cont 34 A very interesting and useful tool for a 3D visualization in GRASS is the nviz tool. At the GRASS prompt you type GRASS 6.3.cvs (w13):~ > nviz v_pred This command will display the following: If you want to display the data points click on "Visualize" and access the data points through "Vector Points". Using the available options you can control different parameters. You can also save this display by clicking on "File" and then on "Save image as" and you can save it as a TIFF file. Exercise: Repeat the above using the data walker lake u.txt. 35 Spatial data analysis with GRASS and R - Exampe We will use the Wolfcamp data from Texas (see Cressie (1993, pp. 212-214). The U.S. Department of Energy proposed (in the 1980s) a nuclear waste site to be in Deaf Smith County in Texas (bordering New Mexico). The contamination of the aquifer was a concern, and therefore the piezometric-head data were obtained at 85 locations by drilling a narrow pipe through the aquifer. The measures are in feet above sea level. See figure below for the location (Texas panhandle): First we need to create a new location that will hold these data. You can name it wf, and also create a mapset a1 . Read the data first in R to find the north-south and east-west extend of the region. To create the location in GRASS we use projection values and choose x, y coordinates. The minimum value of x is -145.2 and its maximum 112.8. The minimum value of y is 9.4 and its maximum 184.8. Therefore define the location as follows: North: 186 South: 8 West: -144 East: 114 and 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 access the ascii file wolfcamp.txt using the v.in.ascii command. Give an output name wolfcamp in. See the snapshots of the v.in.ascii module below: 36 Define the columns of the data set: Once you have imported the file, display the points (first open a monitor window, and then type d.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/circle to get the following plot: 37 From GRASS to R: Now we want to take these data into R, analyze them using gstat, and then get the kriging predictions again back to GRASS. Simply at the command line of GRASS we type R. Once you are in 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 level 1 (42.8, 127.6, 0) 1 42.8 127.6 1464 2 (-27.4, 90.8, 0) 2 -27.4 90.8 2553 3 (-1.2, 84.9, 0) 3 -1.2 84.9 2158 4 (-18.6, 76.5, 0) 4 -18.6 76.5 2455 5 (96.5, 64.6, 0) 5 96.5 64.6 1756 We observe that the data are read as "SpatialPointsDataFrame" (the number in parenthesis are read as coordinates). Now, we can compute variograms, fit variogram models, and do kriging predictions. For this example 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.fit model psill range 1 Nug 9959.434 0.00000 2 Sph 33852.314 66.08075 Here is the fitted model variogram to the sample variogram: 38 We 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 184 grd <- 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+y gridded(grd) <- TRUE class(grd) plot(grd, cex=0.5) points(a, pch=1, col="green", cex=0.9) The grid with the data points: 39 Now 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.var From 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 x0 d.rast level_pred d.rast level_var You 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=50 Reading data: 100% Range of data: min = 785.597107 max = 3647.292236 Range of levels: min = 800.000000 max = 3600.000000 Displacing data: 100% Total levels: 57 WARNING: 97 crossings founds Building topology ... 66 primitives registered Building areas: 100% 0 areas built 0 isles built Attaching islands: Attaching centroids: 100% Topology was built. Number of nodes : 66 Number of primitives: 66 Number of points : 0 Number of lines : 66 Number of boundaries: 0 Number of centroids : 0 Number of areas : 0 Number of isles : 0 40 r.contour complete. Here is the raster map with the contours and the data points: GRASS 6.3.cvs (wf):~ > d.rast level_pred 100% GRASS 6.3.cvs (wf):~ > d.vect level_pred_cont GRASS 6.3.cvs (wf):~ > d.vect wolfcamp_in size=5 icon=basic/circle 41 A very interesting and useful tool for a 3D visualization in GRASS is the nviz tool. At the GRASS prompt you type GRASS 6.3.cvs (wf):~ > nviz level_pred This 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 following three-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 file. 42 Create a LOCATION using a georefrenced map If the map we want to import is georefernced (it has real coordinates - longitude and latitude) we first create a LOCATION by clicking on "Georeferenced file" on the GRASS startup window to open the following dialog box: Give a name to the LOCATION and browse to find the path to the georeferenced map. For this tutorial we are using the map o34118a4.tiff obtained from http://atlas.ca.gov/download.html?sl=casil/imageryBaseMapsLandCover/baseMaps/drg/ The map o34118a4.tiff is a GeoTiff file and it was created using the Albers projection system. 43 We have created a new LOCATION named UCLA WLA as shown on the GRASS startup screen below: We observe that the PERMANENT mapset was created automatically by GRASS, but let's create and enter GRASS through a mapset named ucla1: 44 At the command window we read the file, list the available raster maps (here there is only one - ucla), and open a monitor to display the map: GRASS 6.3.cvs (UCLA_WLA):~ > r.in.gdal in=~/grass/o34118a4.tiff out=ucla Projection of input dataset and current location appear to match. Proceeding with import... 100% GRASS 6.3.cvs (UCLA_WLA):~ > g.list rast ---------------------------------------------<raster> files available in mapset <ucla1>: ucla ---------------------------------------------GRASS 6.3.cvs (UCLA_WLA):~ > d.mon x0 using default visual which is TrueColor ncolors: 16777216 Graphics driver [x0] started GRASS 6.3.cvs (UCLA_WLA):~ > d.rast ucla 100% GRASS 6.3.cvs (UCLA_WLA):~ > Here is the map displayed through the monitor x0: Another way of course to display the map is through the GIS Manager. 45 Let's check the region for this map using the g.region command: GRASS 6.3.cvs (UCLA_WLA):~ > g.region -pec projection: 99 (Albers Equal Area) zone: 0 datum: nad83 ellipsoid: grs80 north: -431016.92261189 south: -445071.36291951 west: 138267.95234148 east: 150034.7091132 nsres: 1.52715857 ewres: 1.52715857 rows: 9203 cols: 7705 cells: 70909115 north-south extent: 14054.440308 east-west extent: 11766.756772 center easting: 144151.330727 center northing: -438044.142766 GRASS 6.3.cvs (UCLA_WLA):~ > We have the following: The projection used for this map is Albers Equal Area, datum is nad83, and ellipsoid is grs80. The numbers north: south: west: east: -431016.92261189 -445071.36291951 138267.95234148 150034.7091132 define the extent of our region. The resolution is 1.52715857 on both n-s and e-s directions. From these we can get the number of rows and columns using the following calculation: Number of rows=(445071.36291951-431016.92261189)/1.52715857=9203 and Number of columns=(150034.7091132-138267.95234148)/1.52715857=7705 The resolution can be changed, as discussed in the previous handout using the command g.region with the options nsres and ewres. For example GRASS 6.3.cvs (UCLA_WLA):~ > g.region nsres=2 ewres=2 46 If we want to see the longitude and latitude for this region we use the flag -l: GRASS 6.3.cvs (UCLA_WLA):~ > g.region -lpec projection: 99 (Albers Equal Area) zone: 0 datum: nad83 ellipsoid: grs80 north: -431016.92261189 south: -445071.36291951 west: 138267.95234148 east: 150034.7091132 nsres: 1.52715857 ewres: 1.52715857 rows: 9203 cols: 7705 cells: 70909115 north-west corner: long: 118:30:03.164016W lat: 34:07:36.185959N north-east corner: long: 118:22:23.973857W lat: 34:07:29.891711N south-east corner: long: 118:22:33.369749W lat: 33:59:53.850963N south-west corner: long: 118:30:11.8234W lat: 34:00:00.135983N center longitude: 118:26:18.082755W center latitude: 34:03:45.016154N north-south extent: 14054.440308 east-west extent: 11766.756772 center easting: 144151.330727 center northing: -438044.142766 GRASS 6.3.cvs (UCLA_WLA):~ > 47 Connecting more than one quadrangles together (by Jan de Leeuw) Again we obtained 6 quadrangles as shown on the map below from http://atlas.ca.gov/casil/imageryBaseMapsLandCover/baseMaps/drg/ 7.5_minute_series_albers_nad83_untrimmed/ http://atlas.ca.gov/casil/imageryBaseMapsLandCover/baseMaps/drg/ 7.5_minute_series_albers_nad83_untrimmed/34118/ We want to import and connect the 6 quadrangles named: Beverly Hills, Van Nuys, Canoga Park, Calabasas, Malibu Beach, and Topanga. The file names of these quadrangles are 034118a4, 034118b4, 034118b5, 034118b6, 034118a6, 034118a5. 48 First let's exit GRASS and enter again through UCLA WLA and the PERMANENT mapset (we will see why...). We then use r.in.gdal to read these files. At the command window we type the following one at a time. r.in.gdal r.in.gdal r.in.gdal r.in.gdal r.in.gdal r.in.gdal in=~/grass/o34118a4.tiff in=~/grass/o34118b4.tiff in=~/grass/o34118b5.tiff in=~/grass/o34118b6.tiff in=~/grass/o34118a6.tiff in=~/grass/o34118a5.tiff out=A out=B out=C out=D out=E out=F And here is what we have... GRASS 6.3.cvs (UCLA_WLA):~ > g.list rast ---------------------------------------------<raster> files available in mapset <PERMANENT>: A B C D E F ---------------------------------------------GRASS 6.3.cvs (UCLA_WLA):~ > We will now create and save the new region that will hold all 6 maps: GRASS 6.3.cvs (UCLA_WLA):~ > g.region -s rast=A,B,C,D,E,F GRASS 6.3.cvs (UCLA_WLA):~ > g.region -p projection: 99 (Albers Equal Area) zone: 0 datum: nad83 ellipsoid: grs80 north: -417152.28930665 south: -445407.37680891 west: 115026.92199555 east: 150034.7091132 nsres: 1.52713693 ewres: 1.52719047 rows: 18502 cols: 22923 cells: 424121346 GRASS 6.3.cvs (UCLA_WLA):~ > 49 We then open a monitor to display the maps. We type the following commands one at a time, and we also use the flag -o which stands for overlay: d.rast d.rast d.rast d.rast d.rast d.rast -o -o -o -o -o -o map=A map=B map=C map=D map=E map=F The 6 maps are displayed one by one in the monitor and finally we get the picture below: We can improve the result by removing the spaces between the quadrangles. Also, we want to create a new map that contains all these 6 quadrangles. This is done with the following commands: 50 GRASS 6.3.cvs (UCLA_WLA):~ > r.series in=A,B,C,D,E,F out=wla method=maximum Percent complete ... 100% GRASS 6.3.cvs (UCLA_WLA):~ > g.list rast ---------------------------------------------<raster> files available in mapset <PERMANENT>: A B C D E F wla ---------------------------------------------GRASS 6.3.cvs (UCLA_WLA):~ > We have used the command r.series with the option method=maximum which takes the maximum of all pixels in the region over all rasters. The input maps are the 6 quadrangles and we named the output map wla. We also want to keep the original colors (from A, B, C, D, E, F). We can choose any one of the 6 maps to be the base color for our newly created map wla. Here is the command: GRASS 6.3.cvs (UCLA_WLA):~ > r.colors map=wla rast=A Color table for [wla] set to A GRASS 6.3.cvs (UCLA_WLA):~ > Note: Other options for method= average: average value count: count of non-NULL cells median: median value mode: most frequently occuring value minimum: lowest value stddev: standard deviation sum: sum of values etc. 51 We display the map as shown below: GRASS 6.3.cvs (UCLA_WLA):~ > d.rast wla 100% GRASS 6.3.cvs (UCLA_WLA):~ > 52 We can also connect the 6 quadrangles using the r.patch command as follows. Again, don't forget the r.colors command. Here is the r.patch command with the flag -z: GRASS 6.3.cvs (UCLA_WLA):~ > r.patch -z in=A,B,C,D,E,F out=wla2 r.patch: percent complete: 100% Creating support files for wla2 GRASS 6.3.cvs (UCLA_WLA):~ > Finally, instead of typing all these commands, we can type only the module name (e.g. r.series, r.patch), which will take us to the dialog box and interactively select the maps to be connected, etc. Here is the snapshot of the r.series, and r.patch dialog boxes: 53 More commands... To rename a raster map we use g.rename rast=old mapname, new mapname. To remove a raster map we use g.remove rast=mapname. To access raster maps in other mapsets other than the current mapset: Suppose you are in mapset ucla1 of the location UCLA WLA but you want to display the map wla from the PERMANENT mapset. You can type: GRASS 6.3.cvs (UCLA_WLA):~ > d.rast wla@PERMANENT 100% GRASS 6.3.cvs (UCLA_WLA):~ > To export a map as a TIFF file you use: GRASS 6.3.cvs (UCLA_WLA):~ > r.out.tiff in=ucla13 out=/Users/nicolas/Desktop/ucla_area 100% r.out.tiff complete. GRASS 6.3.cvs (UCLA_WLA):~ > You can also use the dialog box after you type the command r.out.tiff. Note: GRASS can export raster data to more than 20 different formats using the r.out.gdal module. Some of them are GTiff, PNG, JPEG, GIF, BMP, etc. 54 More on raster maps Creating a new LOCATION using the r.in.gdal command: The command r.in.gdal is another way to create a location. Suppose that LOCATION west1 already exists. We enter GRASS through west1 and mapset aaaaa. We are at the command window as follows: Welcome to GRASS 6.3.cvs (2007) GRASS homepage: http://grass.itc.it/ This version running thru: Bash Shell (/bin/bash) Help is available with the command: g.manual -i See the licence terms with: g.version -c If required, restart the graphical user interface with: gis.m When ready to quit enter: exit GRASS 6.3.cvs (west1):~ > To create a new location at the command line we type: GRASS 6.3.cvs (west1):~ > r.in.gdal ~/grass/snow.jpg > out=snow 100% 100% 100% GRASS 6.3.cvs (west1):~ > location=west1a What happened here? r.in.gdal reads the file snow.jpg from the home directory /grass, gives the output name snow and creates a new location names west1a. The coordinates of the new location west1a will be the same as with the location west1. To check what we just did, exit from GRASS and go back to the startup screen: 55 We see the new location west1a and its PERMANENT mapset. We can enter GRASS through this location and check the coordinates system of west1a and also list the raster maps. GRASS 6.3.cvs (west1a):~ > g.region -pec projection: 0 (x,y) zone: 0 north: 1878 south: 0 west: 0 east: 900 nsres: 1 ewres: 1 rows: 1878 cols: 900 cells: 1690200 north-south extent: 1878.000000 east-west extent: 900.000000 center easting: 450.000000 center northing: 939.000000 GRASS 6.3.cvs (west1a):~ > g.list rast ---------------------------------------------<raster> files available in mapset <PERMANENT>: snow.blue snow.green snow.red ---------------------------------------------GRASS 6.3.cvs (west1a):~ > More than one monitors can be opened from the command window. GRASS can open up to 7 windows (x0, x1, x2, x3, x4, x5, x6). You can open them as follows: GRASS 6.3.cvs (west1a):~ > d.mon x0 using default visual which is TrueColor ncolors: 16777216 Graphics driver [x0] started GRASS 6.3.cvs (west1a):~ > d.mon x1 using default visual which is TrueColor ncolors: 16777216 Graphics driver [x1] started GRASS 6.3.cvs (west1a):~ > d.mon x2 using default visual which is TrueColor ncolors: 16777216 Graphics driver [x2] started GRASS 6.3.cvs (west1a):~ > The current monitor is the last one opened (here x2). To switch to another monitor (e.g. to x0) you can type: GRASS 6.3.cvs (west1a):~ > d.mon select=x0 56 ...
View Full Document

This note was uploaded on 02/11/2012 for the course STATS c173/c273 taught by Professor Nicolaschristou during the Spring '11 term at UCLA.

Ask a homework question - tutors are online