Difference between revisions of "Software"

m (Replaced content with "{{#meta: REFRESH | 1;url=http://gsif.isric.org/doku.php/wiki:software}} #REDIRECTS to [http://gsif.isric.org/doku.php/wiki:software GSIF tutorials]")
Line 1: Line 1:
{{Drop|T}}his article provides some basic instructions on how to install and do first steps in software the R + Open Source Desktop GIS (OSGeo) + Google Earth bundle that can be used to run analysis of various spatio-temporal data (this article was specifically written for the MS Windows OS users). This bundle allows full GIS+statistics integration and can support practically about 80% of processing/visualization capabilities available in commercial packages such as ArcInfo/Map or Idrisi (read more about the [http://r4stats.com/articles/popularity/ popularity of R]). Although R contains an increasing number of packages for processing GIS raster layers e.g. [http://cran.r-project.org/web/packages/biOps/ biOps], [http://cran.r-project.org/web/packages/rimage/ rimage]; [http://r-forge.r-project.org/projects/raster/ raster] (see also the list of [http://r-forge.r-project.org/projects/r-gis r-gis] projects under development), the real power is in combining R with Open Source GIS such as SAGA GIS and GRASS GIS.
{{#meta: REFRESH | 1;url=http://gsif.isric.org/doku.php/wiki:software}}
#REDIRECTS to [http://gsif.isric.org/doku.php/wiki:software GSIF tutorials]
Software such as SAGA GIS allow you to input, edit and visually explore geographical data, before and after the actual statistical analysis. Note that there is a large number of [http://freegis.org/ open GIS] alternatives and commercial GIS packages that could also be linked to R with equal success. For example, GRASS GIS (versions [http://grass.osgeo.org/ >6.3] also available for MS Windows machines) is another good candidate for GIS+statistics integration. Within GRASS, you can load R and then use its functionality on the GRASS map layers, or run R on top of GRASS (with the help of the [http://cran.r-project.org/web/packages/spgrass6/ spgrass6] package). Another alternative for GIS+R integration is the [http://qgis.org/ QGIS], which has among its main characteristics a python console, and a very elaborated way of adding Python plugins, which is already used for an efficient R plugin ([http://www.ftools.ca manageR]).
[[Image:R_SAGA_GE_logo.jpg|right|thumb|365px|The R + Open Source desktop GIS + Google Earth bundle scheme.]]
= R =
== Installation ==
To install R under Windows, download and run an '''[http://cran.r-project.org/ installation exe]''' from the R project homepage. This will install R for window with a GUI. After you start R, you will first need to set-up the working directory and install additional packages. To run geostatistical analysis in R, you will need to add the following R packages: [http://www.gstat.org/ gstat] (gtat in R), [http://rgdal.sourceforge.net/ rgdal] (GDAL import of GIS layers in R), [http://r-spatial.sourceforge.net/ sp] (operations on maps in R), foreign, [http://www.spatstat.org/ spatstat] (spatial statistics in R) and [http://cran.r-project.org/web/packages/maptools/index.html maptools] (see also: a complete list of [http://cran.r-project.org/web/packages/ '''contributed packages''']). To install these packages you should do the following. First start the R GUI, then select ''Packages > Load package'' from the main menu. Note that, if you wish to install a package on the fly, you will need to select a suitable CRAN mirror from where it will download and unpack a package. You can also install a number of packages used for the [http://cran.r-project.org/web/views/Spatial.html spatial analysis] by typing e.g.:
<geshi lang=R lines=0>
> install.packages(c("rgdal","maptools","gstat","spdep","spatstat","RSAGA","spgrass6"))
Another quick way to get all packages used in R to do spatial analysis (as explained in the [http://www.asdar-book.org/ ASDAR book] by [http://www.springer.com/978-0-387-78170-9 Bivand et al. 2008]) is to install the [http://cran.r-project.org/web/packages/ctv/index.html ctv] package and then execute the command:
<geshi lang=R lines=0>
> install.packages("ctv")
> library(ctv)
> install.views("Spatial")
If you have problems accessing the web, you can also try downloading all packages of interest from CRAN (see e.g.: [http://cran.r-project.org/bin/windows/contrib/r-release/ windows '''precompiled binaries''']) and then install them from the local file. After you install a package, you will still need to load it into your workspace (every time you start R) before you can use its functionality. A package is commonly loaded using e.g.:
<geshi lang=R lines=0>
> library(gstat)
== First steps in R ==
R is today identified as one of the [http://r4stats.com/articles/popularity/ fastest growing] and most comprehensive statistical computing tools/communities. It practically offers statistical analysis and visualisation of unlimited sophistication. A user is not restricted to a small set of procedures or options, and because of the contributed packages, users are not limited to one method of accomplishing a given computation or graphical presentation ([http://www.springer.com/978-0-387-78170-9 Bivand et al., 2008]; [http://www.itc.nl/~rossiter/teach/R/RIntro_ITC.pdf Rossiter, 2007]; [http://www.crcpress.com/shopping_cart/products/product_detail.asp?sku=C486X Murrell, 2006]). R became attractive for analysis of spatial data mainly due to the recent integration of the geostatistical tools (e.g. [http://www.gstat.org gstat], [http://leg.ufpr.br/geoR/ geoR]) and tools that allow R computations with spatial data layers (e.g. [http://r-spatial.sourceforge.net/ sp/maptools], [http://rgdal.sourceforge.net/ rgdal], [http://r-gis.r-forge.r-project.org/ r-gis]). A complete [http://cran.r-project.org/web/views/Spatial.html listing] of all R packages that can be used to run analysis on spatial data in R is maintained by Roger Bivand. The '[http://cran.r-project.org/web/views/Spatial.html Spatial]' packages can be nicely combined with e.g. the [http://cran.r-project.org/web/views/Environmetrics.html Environmentrics] packages. The [http://www.interactivegraphics.org interactive graphics] in R is also increasingly powerful.
Although there are few GUI's being developed to run analysis in R (e.g. [http://socserv.mcmaster.ca/jfox/Misc/Rcmdr/ Rcmdr]), R is (and will stay) a command-based software. So, the first step you should consider is to get some initial training on [[Scripting in R|R syntax]]. You may consult some literature first, but also try to follow 1-2 days training courses. The second useful step to get into R is to try to obtain a more user-friendly R scripting edition such as as '''[http://www.sciviews.org/Tinn-R/ TINN-R]''', '''[http://www.rstudio.org/ Rstudio] and/or '''[http://rosuda.org/JGR/ JGR]'''. You can use this [http://spatial-analyst.net/scripts/Rprofile.site Rprofile.site] (copy to your "R/etc/" directory) and this [http://spatial-analyst.net/scripts/Rconsole Rconsole] setting to produce the same GUI configuration as shown below.
If a single argument in the command is incorrect, inappropriate or mistyped, you will get an error message. If the error message is not helpful, then try receiving more help about some operation by typing e.g.:
<geshi lang=R lines=0>
> help(krige)
<geshi lang=R lines=0>
> ?krige
commands or via the Html help files. You may also wish to register to the special interest groups such as [https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo] or similar and subscribe to their mailing lists. [http://mailman.geo.uu.nl/mailman/listinfo/gstat-info Gstat-info] is the mailing list of the [http://www.gstat.org/ gstat] package and also offers many interesting information.
These are some easy to follow R tutorials:
# [http://www.math.csi.cuny.edu/Statistics/R/simpleR/ Simple R]
# [http://www.statmethods.net/ Quick R]
# [http://www.fort.usgs.gov/brdscience/learnRE.htm Short course in R] by Paul Geissler (USGS)
# [http://www.itc.nl/personal/rossiter/teach/R/RIntro_ITC.pdf Introduction to R] by D.G. Rossiter
# [http://www.bias-project.org.uk/ASDARcourse/ Introductory course on Spatial Data Analysis with R] by R. Bivand
A brief guide to learn how to make scripts in R is available in this [[Scripting in R|article]]. But please, consider also purchasing some of the official [http://www.r-project.org/doc/bib/R-books.html R books] and/or attending some of the short [[Training in R|courses on R]].
To download the file ([http://geomorphometry.org/content/baranja-hill DEM25m]) used in this exercise, we can do:
<geshi lang=R lines=0>
> download.file("http://spatial-analyst.net/DATA/DEM25m.zip", destfile=paste(getwd(),"DEM25m.zip",sep="/"))
> fname <- zip.file.extract(file="DEM25m.asc", zipname="DEM25m.zip")
> file.copy(fname, "./DEM25m.asc", overwrite=TRUE)
Next, we can import the ArcInfo ASCII map to R using the rgdal package:
<geshi lang=R lines=0>
> library(rgdal)
> dem25m <- readGDAL("DEM25m.asc")
After we attach the correct coordinate system string to this map, we can export it to e.g. ILWIS GIS by using:
<geshi lang=R lines=0>
> proj4string(dem25m) <- CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel
  +    +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m")
> writeGDAL(dem25m[1], "dem25m.mpr", "ILWIS")
A serious limitation of R, considering GIS data, is the memory handling issue. You will soon discover that it is almost impossible to visualize and process raster maps which are >1--5 million pixels (see also the [http://raster.r-forge.r-project.org/ raster project]). In addition, R has still limited capabilities for visual data exploration, and almost no capabilities for data input/editing. You can overcome these problems by combining R with the open source GIS, such as [[Software#SAGA GIS|SAGA]], [http://www.osgeo.org/files/journal/final_pdfs/OSGeo_vol1_GRASS-R.pdf GRASS] and/or [[Software#ILWIS GIS|ILWIS]].
== Installing SAGA ==
[[Image:Fig_R_SAGA_screen.jpg|thumb|450px|Execution of SAGA command line via R and example of an output window.]]
SAGA stands for [http://saga-gis.org '''System for Automated Geoscientific Analyses'''] is a full-fledged GIS, and many of its features have some relation with geomorphometry, which makes it an ideal tool for operational work, but also for GIS training purposes. It is an open-source GIS with support for raster and vector data. It includes a large set of geoscientific algorithms, and is especially powerful for the analysis of DEMs.
SAGA has been under development since 2001 at the University of Göttingen (the SAGA development team, including the key developer Olaf Conrad, has since moved to University of Hamburg), Germany, with aim of simplifying the implementation of new algorithms for spatial data analysis within a framework that immediately allows their operational application. Therefore, SAGA targets not only the pure user but also the developer of geo-scientific methods. SAGA has its roots in DiGeM, a small program specially designed for the extraction of hydrological land-surface parameters. In 2004 most of SAGA's source code was published using an Open Source Software (OSS) license. To learn more about SAGA, try obtaining some of the following [http://www.saga-gis.org/en/about/references.html references]:
*Böhner, J., McCloy, K.R., Strobl, J. (Eds.), 2006. [http://www.saga-gis.org/en/about/references.html '''SAGA — Analysis and Modelling Applications''']. Göttinger Geographische Abhandlungen, Heft 115. Verlag Erich Goltze GmbH, Göttingen, 117 pp.
*Brenning, A. 2008. [http://downloads.sourceforge.net/saga-gis/hbpl19_03.pdf Statistical geocomputing combining R and SAGA: The example of landslide susceptibility analysis with generalized additive models]. In: J. Böhner, T. Blaschke & L. Montanarella (eds.), SAGA - Seconds Out (Hamburger Beiträge zur Physischen Geographie und Landschaftsökologie, 19), 23-32.
*Conrad, O. 2007. [http://webdoc.sub.gwdg.de/diss/2007/conrad/conrad.pdf '''SAGA - Entwurf, Funktionsumfang und Anwendung eines Systems für Automatisierte Geowissenschaftliche Analysen''']. electronic doctoral dissertation, University of Göttingen.
*Olaya, V., Conrad, O., 2008. [http://dx.doi.org/10.1016/S0166-2481(08)00012-3 '''Geomorphometry in SAGA''']. In: Hengl, T. and Reuter, H.I. (Eds), Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, vol. 33, Elsevier, 293-308 pp.
SAGA can be obtained from a Source Forge [http://sourceforge.net/projects/saga-gis/ repository]; the most recent version is available both for windows and linux operating systems. The fastest installation is that you simply obtained the compiled binaries (e.g. 'saga_2.*.*_*.zip') and unzip them into the "/library/RSAGA/saga_vc/" or to "C:/saga_vc" because this is the default location where RSAGA looks for the binary files. Here you can obtain the list of all [http://sourceforge.net/projects/saga-gis/ SAGA GIS (2.1.2)] modules available via the SAGA command line. For a complete description and usage of function arguments see the [http://spatial-analyst.net/DATA/SAGA_GIS_2_1_2.txt text file] or refer to the [http://sourceforge.net/apps/trac/saga-gis/ SAGA help documentation].
Important sources:
* [http://spatial-analyst.net/DATA/SAGA_GIS_2_1_2.txt Complete list of SAGA GIS modules] (with parameters)
* [http://sourceforge.net/apps/trac/saga-gis/ The official SAGA help] documentation page
== First steps in SAGA ==
First, you need to load the [http://cran.r-project.org/web/packages/RSAGA/index.html RSAGA] library by:
<geshi lang=R lines=0>
> library(RSAGA)
Note that loading a library you will not be able to run any of the SAGA operations if also SAGA binaries and libraries are not installed on the same machine (SAGA is an external GIS packages and NOT an built-in R package). You can at any time set the location of your SAGA environment by typing:
<geshi lang=R lines=0>
> myenv <- rsaga.env(path="C:/saga_vc")
To find out what a certain module (the best way to browse all possible modules is to open SAGA and look at the left menu tab that shows all [http://www.saga-gis.org/saga_modules_doc/index.html module libraries] available) in SAGA does, you can execute:
<geshi lang=R lines=0>
> rsaga.get.modules("grid_gridding")
which will show this:
<geshi lang=R lines=0>
#! $grid_gridding
#!  code                      name interactive
#! 1    0            Shapes to Grid      FALSE
#! 2    1 Inverse Distance Weighted      FALSE
#! 3    2        Nearest Neighbour      FALSE
#! 4    3        Natural Neighbour      FALSE
#! 5    4 Modifed Quadratic Shepard      FALSE
#! 6    5            Triangulation      FALSE
#! 7    6 Kernel Density Estimation      FALSE
#! 8    7 Angular Distance Weighted      FALSE
#! Warning message:
#! NAs introduced by coercion
To get further usage of the 3rd function of the grid_gridding module, you can execute:
<geshi lang=R lines=0>
> rsaga.get.usage("grid_gridding", 2)
which then gives the complete list of parameters (obligatory and optional) that need to be set before this command can be executed:
<geshi lang=R lines=0>
#! library path:  C:\saga_vc\modules\grid_gridding.dll
#! library name:  Grid - Gridding
#! Usage: saga_cmd grid_gridding 0 -INPUT <str> [-FIELD <str>] [-OUTPUT <str>] [-MULTIPLE <str>] [-LINE_TYPE <str>] [-POLY_TYPE <str>]
#!  [-GRID_TYPE <str>] [-TARGET <str>] [-USER_XMIN  <str>] [-USER_XMAX <str>] [-USER_YMIN <str>] [-USER_YMAX <str>] [-USER_SIZE <str>]
#!  [-USER_FIT <str>] [-USER_GRID <str>] [-USER_COUNT <str>] [-USER_BCOUNT] [-GRID_GRID <str>] [-GRID_COUNT <str>]
#!  -INPUT:<str>          Shapes
#!        Shapes (input)
#!  -FIELD:<str>          Attribute
#!        Table field
#!  -OUTPUT:<str>        Output Values
#!        Choice
#!        Available Choices:
#!        [0] data / no-data
#!        [1] index number
#!        [2] attribute
#!        Default: 2
Another place to look for libraries and parameters is the [http://www.saga-gis.org/saga_modules_doc/ SAGA's homepage] or this [[list of SAGA modules]]. SAGA also provides an [http://www.saga-gis.org/saga_api_doc/html/ API functionality], so you can access SAGA grids and libraries from e.g. python.
Once you have learned about a specific function, you can run it via the ''rsaga.geoprocessor''. For example to extract a stream network from a DEM, we can run:
<geshi lang=R lines=0>
> rsaga.esri.to.sgrd(in.grids="dem25m.asc", out.sgrds="dem25m.sgrd", in.path=getwd())
> rsaga.geoprocessor(lib="ta_channels", module=0, param=list(ELEVATION="dem25m_f.sgrd",
  +    CHNLNTWRK="channel_ntwrk.sgrd", CHNLROUTE="channel_route.sgrd", SHAPES="channels.shp",
  +    INIT_GRID="dem25m_f.sgrd", DIV_CELLS=3, MINLEN=40))
> channels <- readOGR("channels.shp", "channels")
> spplot(channels["LENGTH"], col.regions=bpy.colors())
which will produce a screen as shown on the right (above). Note that the ''rsaga.geoprocessor'' only serves to pass the parameters to SAGA command line (saga_cmd.exe). If you are loading a list of grids to some SAGA process, they need to be written as a list:
<geshi lang=R lines=0>
# export grids back to ArcInfo ASCII format:
> rsaga.sgrd.to.esri(in.sgrds=c("channel_ntwrk.sgrd","channel_route.sgrd"),
+      out.grids=c("channel_ntwrk.asc","channel_route.asc"), out.path=getwd(), prec=1)
the rsaga.geoprocessor would expect grids separated by ";", e.g. GRIDS="DSM1.sgrd;DSM2.sgrd;DSM3.sgrd". Sometimes you might receive an error message or computation will fail, which possibly means a bug in SAGA GIS. Such problems you should report directly via the [http://sourceforge.net/forum/forum.php?forum_id=790705 SAGA user's forum] (and not via the R-sig-geo mailing list). Read [[Analysis of DEMs in R+ILWIS/SAGA|'''more''']] on how to extract land-surface parameters and automate analysis of DEMs using SAGA.
== Installing GRASS (MS Windows) ==
[[Image:R_GRASS_screen.jpg|thumb|450px|Extraction of stream network using GRASS functionality from R.]]
GRASS GIS, now one of the eight initial Software Projects of the Open Source Geospatial Foundation, is probably the most known open source GIS software and does need to be specially introduced. Its functionality and usage are described in detail in [http://www.grassbook.org Neteler et al. (2008)]; a wiki-based tutorial is available at http://grass.osgeo.org/wiki/. GRASS itself is a collection of environment variables (they vary from version to version), which are populated when GRASS is started --- that is a GRASS shell is started within the running command shell. GRASS is a much larger projected than SAGA considering the number of developers/institutions involved, although their functionality considering the DEM analysis is about similar. The latest version of GRASS (6.3.) contains over 350 routines.
Although originally a Linux-based project, the most recent version of GRASS (6.3; development version) is now also available for MS Windows machines. The coming version 6.4. is expected to be even more stable under Windows. To install GRASS, simply obtain the [http://grass.osgeo.org/download/ installation exe] and install GRASS to some default location e.g. "C:/GRASS/"
== Controlling GRASS from R ==
The GRASS+R integration requires the [http://cran.r-project.org/web/packages/spgrass6/ spgrass6] package. This needs to be installed with all necessary dependencies:
<geshi lang=R lines=0>
> install.packages("spgrass6", dependencies=T)
> library(spgrass6)
Next, you need to set the location of your GRASS installation:
<geshi lang=R lines=0>
> loc <- initGRASS("C:/GRASS", home=tempdir())
> loc
#! gisdbase    c:/Temp
#! location    file678418be
#! mapset      file3d6c4ae1
#! rows        1
#! columns    1
#! north      1
#! south      0
#! west        0
#! east        1
#! nsres      1
#! ewres      1
#! projection  NA
This will automatically generate a temporary GRASS gisdbase (and create a temporary file). Note that this database contains no maps. These can be imported using:
<geshi lang=R lines=0>
> execGRASS("r.in.gdal", flags="o", parameters=list(input="DEM25m.asc", output="DEM"))
#! WARNING: Over-riding projection check
#! RINGDA~1 complete.
> execGRASS("g.region", parameters=list(rast="DEM"))
> gmeta6()
#! gisdbase    c:/Temp
#! location    file678418be
#! mapset      file3d6c4ae1
#! rows        149
#! columns    147
#! north      5074275
#! south      5070550
#! west        6551872
#! east        6555547
#! nsres      25
#! ewres      25
#! projection  NA
To get a description of some GRASS function and parameters connected, you can use e.g.:
<geshi lang=R lines=0>
> parseGRASS("r.in.gdal")
Once you have imported a map to GRASS format, you can also open and visualize this DEM using the GRASS GUI. In the next step, we can extract the drainage network:
<geshi lang=R lines=0>
> execGRASS("r.watershed", flags="overwrite",
+    parameters=list(elevation="DEM", stream="stream", threshold=as.integer(50)))
#! SECTION 1a (of 5): Initiating Memory.
#! SECTION 1b (of 5): Determining Offmap Flow.
#! SECTION 2: A * Search.
#! SECTION 3: Accumulating Surface Flow.
#! SECTION 4: Watershed determination.
#! SECTION 5: Closing Maps.
# thin the raster map so it can be converted to vectors:
> execGRASS("r.thin", parameters=list(input="stream", output="streamt"))
#! File stream -- 149 rows X 147 columns
#! Bounding box: l = 2, r = 148, t = 2, b = 150
#! Pass number 1
#! Deleted 105 pixels
#! Pass number 2
#! Deleted 0 pixels
#! Thinning completed successfully.
#! Output file 149 rows X 147 columns
#! Window 149 rows X 147 columns
# convert to vectors:
> execGRASS("r.to.vect", parameters=list(input="streamt", output="streamt", feature="line"))
#! WARNING: Default driver / database set to:
#!          driver: dbf
#!          database: $GISDBASE/$LOCATION_NAME/$MAPSET/dbf/
#! Extracting lines...
#! Building topology ...
#! Registering lines:
#! 287 primitives registered
#! Building areas:
#! 0 areas built
#! 0 isles built
#! Attaching islands:
#! Attaching centroids:
#! Topology was built.
#! Number of nodes    :  266
#! Number of primitives:  287
#! Number of points    :  0
#! Number of lines    :  287
#! Number of boundaries:  0
#! Number of centroids :  0
#! Number of areas    :  0
#! Number of isles    :  0
#! RTOVEC~1 complete.
Finally, we can import the GRASS maps to R using the spgrass6 package:
<geshi lang=R lines=0>
> streamt <- readVECT6("streamt")
#! Exporting 287 points/lines...
#! 287 features written
#! OGR data source with driver: ESRI Shapefile
#! Source: "c:/Temp/file678418be/file3d6c4ae1/.tmp", layer: "streamt"
#! with  287  rows and  3  columns
#! Feature type: wkbLineString with 2 dimensions
> plot(streamt)
= FWTools =
[http://fwtools.maptools.org FWTools] (created and maintained by [http://home.gdal.org/warmerda/ Frank Warmerdam]; director and active contributor to OSGeo) is a stand-alone collection of open source tools for processing spatial data. It includes OpenEV, GDAL, MapServer, PROJ.4 and OGDI --- i.e. it contain a number of (highly efficient!) [http://www.gdal.org/gdal_utilities.html utilities] excellent for processing large datasets. Or to quote the author: FWTools ''is intended to give folks a chance to use the latest and greatest functionality''.
== Installation ==
Linux or Windows version of FWTools can be obtained from the [http://fwtools.maptools.org maptools.org] website. After you install the software, you can test it by starting the OpenEV and opening and overlaying some GIS layers. FWTools does not much GUI --- we are mainly interested in using the [http://www.gdal.org/gdal_utilities.html GDAL utilities]. The complete description of the package functionality is available via the [http://trac.osgeo.org/gdal/wiki/FWTools GDAL Wiki]. There is also the [http://lists.maptools.org/mailman/listinfo/fwtools FWTools mailing list].
== Running FWTools from R  ==
There is still no package to control FWTools from R, but we can simply send command lines using the <tt>system</tt> command. Before we can use [http://fwtools.maptools.org FWTools] from R, we need to locate it on our PC by using e.g.:
<geshi lang=R lines=0>
> fw.path = utils::readRegistry("SOFTWARE\\WOW6432Node\\FWTools")$Install_Dir
> gdalwarp = shQuote(shortPathName(normalizePath(file.path(fw.path, "bin/gdalwarp.exe"))))
> gdalwarp
#! [1] "C:\\PROGRA~2\\FWTOOL~1.7\\bin\\gdalwarp.exe"
> workd <- paste(gsub("/", "\\\\", getwd()), "\\", sep="")
Now we can download some GIS data from web:
<geshi lang=R lines=0>
> MOD12Q1 <- "ftp://anonymous:test@e4ftl01u.ecs.nasa.gov/
+        MOLT/MOD12Q1.004/2004.01.01/"
> download.file(paste(MOD12Q1, "MOD12Q1.A2004001.h18v03.004.2006117173748004.2006117173748.hdf", sep=""),
+        destfile=paste(getwd(),
+        "MOD12Q1.A2004001.h18v03.004.2006117173748.hdf", sep="/"), mode='wb')
#! Resolving e4ftl01u.ecs.nasa.gov...
#! Connecting to e4ftl01u.ecs.nasa.gov||:21... connected.
#! Logging in as anonymous ... Logged in!
#! ==> SYST ... done.    ==&gt; PWD ... done.
#! ==> TYPE I ... done.  ==&gt; CWD /MOLT/MOD12Q1.004/2004.01.01 ... done.
#! ==> SIZE MOD12Q1.A2004001.h18v03.004.2006117173748.hdf ... 23165983
#! ==> PASV ... done.    ==&gt; RETR MOD12Q1.A2004... done.
#! Length: 23165983 (22M)
#!      0K ..........  ..........  0% 64.9K 5m48s
#!    ...
#!  22550K ..........  .......... 99%  501K 0s
#!  22600K ..........            100%  503K=65s
We can reproject/resample the map to our local coordinate system using the [http://www.gdal.org/gdalwarp.html gdalwarp] functionality (this combines several processing steps in one function):
<geshi lang=R lines=0>
# Define the target coordinate system:
> NL.prj <- "+proj=sterea +lat_0=52.15616055555555
+        +lon_0=5.38763888888889 +k=0.999908 +x_0=155000
+        +y_0=463000 +ellps=bessel +units=m +no_defs
+        +towgs84=565.237,50.0087,465.658,
+          -0.406857,0.350733,-1.87035,4.0812"
# reproject and resample the MODIS image:
> system(paste(gdalwarp, " HDF4_EOS:EOS_GRID:\"", workd,
+        "\\MOD12Q1.A2004001.h18v03.004.2006117173748.hdf\":MOD12Q1:Land_Cover_Type_1
+        -t_srs \"", NL.prj, "\" IGBP2004NL.tif -r near -te 0 300000 280000 625000 -tr 500 500", sep=""))
#! Creating output file that is 560P x 650L.
#! Processing input file HDF4_EOS:EOS_GRID:\\MOD12Q1.A2004001...
#! Using internal nodata values (eg. 255) for image HDF4_EOS:EOS_...
#! 0...10...20...30...40...50...60...70...80...90...100 - done.
In this case we have produced a [https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/land_cover/yearly_l3_global_0_05deg_cmg/mod12c1 MODIS-based land cover map] for the whole Netherlands in resolution of 500 m (in local coordinate system).
In summary, GDAL utilities and similar functions available via FWTools are extremely compact and efficient, and fit to work with large datasets. Commands such as [http://www.gdal.org/gdalwarp.html gdalwarp] allow you to combine several GIS operations in one line and hence easily transfer GIS data from one to other format/grid system. The only thing that might slow you down is the time needed to get the commands correct, which can take more than few iterations.

Latest revision as of 07:54, 13 April 2017

  1. REDIRECTS to GSIF tutorials