Просмотр исходного текста страницы Анализ данных с использованием GRASS GIS и R

==Чтение и запись данных GRASS в среде R==

===Первые шаги===

Запустим GRASS и из командной строки GRASS загружаем R. В конечном итоге перед нами окажется командная строка системы R.

Загрузим в R необходимые библиотеки:

 library(spgrass6)

После загрузки данной библиотеки уже можно взаимодействовать с GRASS GIS. Например, введем следующую команду:

 str(gmeta6())

Результат выполнения команды приводится ниже:

 
 &gt; str(gmeta6())
 List of 24
  $ GISDBASE     : chr "/home/dima/GIS/grass"
  $ LOCATION_NAME: chr "grass"
  $ MAPSET       : chr "PERMANENT"
  $ MONITOR      : chr "x0"
  $ GRASS_GUI    : chr "text"
  $ n            : num 57.3
  $ s            : num 49
  $ w            : num 75
  $ e            : num 90
  $ t            : num 1
  $ b            : num 0
  $ nsres        : num 0.00417
  $ nsres3       : num 1.04
  $ ewres        : num 0.00417
  $ ewres3       : num 1
  $ tbres        : num 1
  $ rows         : int 1992
  $ rows3        : int 8
  $ cols         : int 3600
  $ cols3        : int 15
  $ depths       : int 1
  $ cells        : chr "7171200"
  $ 3dcells      : chr "120"
  $ proj4        : chr "+proj=longlat +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000"
  - attr(*, "class")= chr "gmeta6"

Как легко заменить, данная команда вывела информацию о текущей области (region) GRASS и другие параметры.

Находясь в R, вы имеете возможность запускать команды GRASS, не выходя из сессии R. Для этой цели можно использовать команду system, например:

 
 &gt; system("g.region -p")
 projection: 3 (Latitude-Longitude)
 zone:       0
 datum:      wgs84
 ellipsoid:  wgs84
 north:      57:18N
 south:      49N
 west:       75E
 east:       90E
 nsres:      0:00:15
 ewres:      0:00:15
 rows:       1992
 cols:       3600
 cells:      7171200

===Импорт векторных данных===

Импортируем несколько карт из GRASS GIS в среду R, для этого используется команда readVECT6.

 admin &lt;- readVECT6("admin", ignore.stderr=TRUE)

Просмотрим информацию:

 
 &gt; summary(admin)
 Object of class SpatialPolygonsDataFrame
 Coordinates:
        min      max
 x 75.08194 89.86998
 y 49.08308 57.25027
 Is projected: FALSE 
 proj4string :
 [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]
 Data attributes:
       cat                          NAME  
  Min.   :1.00   Алтайский край       :1  
  1st Qu.:1.75   Новосибирская область:1  
  Median :2.50   Республика Алтай     :2  
  Mean   :2.50                            
  3rd Qu.:3.25                            
  Max.   :4.00

Аналогично, загрузим данные по почвам:

 
 &gt; soils &lt;- readVECT6("soils", ignore.stderr=TRUE) 
 &gt; summary(soils)
 Object of class SpatialPolygonsDataFrame
 Coordinates:
        min      max
 x 75.08194 89.86998
 y 49.08308 57.25027
 Is projected: FALSE 
 proj4string :
 [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]
 Data attributes:
       cat         FAOSOIL         UNIT                   TYPE           NAME   
  Min.   :  1   Ch1-3a : 16   I      :39   LITHOSOLS        :39   Черноземы:52  
  1st Qu.: 44   C1-3a  : 14   Ch     :17   Haplic Chernozems:17   Литосоли :39  
  Median : 87   Sm13-3a: 14   C      :16   CHERNOZEMS       :16   Солонцы  :18  
  Mean   : 87   Oe1-a  : 10   Sm     :15   Mollic Solonetz  :15   Сероземы :11  
  3rd Qu.:130   Mo1-2ab:  9   Ck     :10   Calcic Chernozems:10   Гистосоли:10  
  Max.   :173   (Other):110   (Other):76   (Other)          :76   (Other)  :43  
  NA's   :  1   NA's   :  1   NA's   : 1   NA's             : 1   NA's     : 1  
    FAOSOIL_L    
  Min.   : 1.00  
  1st Qu.: 8.00  
  Median :21.00  
  Mean   :20.26  
  3rd Qu.:34.00  
  Max.   :40.00  
  NA's   : 1.00

При необходимости можно отобразить данные на карте:

 plot(admin, axes=TRUE)

Добавим еще один слой:

 plot(soils, add=TRUE, col="red")

===Импорт растровых данных===

Импорт растровых данных производится аналогично. Для этого предназначена команда readRAST6:

 readRAST6(vname, cat=NULL, ignore.stderr = FALSE, NODATA=NULL, plugin=NULL, mapset=NULL, useGDAL=TRUE, close_OK=TRUE)

В первую очередь, нас будет интересовать параметр vname - название растровой карты GRASS. Считаем, к примеру, растр под названием modis:

 
 &gt; modis &lt;- readRAST6("modis", ignore.stderr=TRUE)
 &gt; summary(modis)
 Object of class SpatialGridDataFrame
 Coordinates:
   min  max
 x  75 90.0
 y  49 57.3
 Is projected: FALSE 
 proj4string :
 [+proj=longlat +a=6378137 +rf=298.257223563 +no_defs
 +towgs84=0.000,0.000,0.000]
 Number of points: 2
 Grid attributes:
   cellcentre.offset    cellsize cells.dim
 x          75.00208 0.004166667      3600
 y          49.00208 0.004166667      1992
 Data attributes:
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
       0    4328    5418    7028    7499   32770  209650 
 

===Экспорт данных в GRASS GIS===

Для экспорта растровых данных существует команда:

 writeRAST6(x, vname, zcol = 1, NODATA=NULL, ignore.stderr = FALSE, useGDAL=TRUE, overwrite=FALSE, flags=NULL)

Здесь x - фрейм пространственных данных, который будет сохранен как растровая карта GRASS под именем vname (подробнее см. в справке R).

Аналогично, для работы с векторными данными используется команда:

 writeVECT6(SDF, vname,  v.in.ogr_flags=NULL, ignore.stderr = FALSE)

Здесь SDF - фрейм пространственных данных, который будет сохранен как векторная карта GRASS под именем vname (подробенее см. в справке R).