Embedded Scribd iPaper - Requires Javascript and Flash Player
4 Auxiliary data sources
As mention previously in §2.10, geostatistical techniques increasingly rely on the availability of auxiliary data sources, i.e. maps and images of surface and sub-surface features. This chapter reviews some of the widely known sources of remote sensing and digital cartographic data that are of interest for various geostatistical mapping applications. I will first focus on the freely available global data sets and remotely sensed data, and then give an example of how to download and import to GIS various MODIS products using scripts. A small repository of global maps at 0.1 arcdegree (10 km) resolution is also available from the authors website1 (see further chapter 7).
R
4.1 Global data sets
Global maps of our environment (both remote sensing-based and thematic) are nowadays increasingly attractive for environmental modeling at global and continental scales. A variety of publicly available maps can be obtained at no cost at resolutions up to 1 km or better, so that the issue for mapping teams is not any more whether to use this data, but where to obtain it, how to load it to an existing GIS and where to find necessary metadata. The most well known global data sets (sorted thematically) are: Height/geomorphology data — Global SRTM Digital Elevation Model is possibly the most known global environmental data set (Rabus et al., 2003). The area covered is between 60° North and 58° South. It was recorded by X-Band Radar (NASA and MIL, covering 100% of the total global area) and CBand Radar (DLR and ASI, covering 40%). The non-public DLR-ASI data is available with a resolution of approximately 30 m (1 arcsec). A complete land surface model ETOPO1 Global Relief Model2 (Fig. 4.1; includes bathymetry data) is available at resolution of 1 km and can be obtained from the NOAA’s National Geophysical Data Center (Amante and Eakins, 2008). The 90 m SRTM DEMs can be obtained from the CGIAR3 — Consortium for Spatial Information. An updated 1 km resolution global topography map (SRTM30 PLUS4 ; used by Google Earth) has been prepared by Becker et al. (2009). From June 2009, ASTER-based Global Digital Elevation Model (GDEM) at resolution of 30 m has been made publicly available. The GDEM was created by stereo-correlating the 1.3 million-scene ASTER archive of optical images, covering almost 98% of Earth’s land surface (Hayakawa et al., 2008). The one-by-one-degree tiles can be downloaded from NASA’s EOS data archive and/or Japan’s Ground Data System5 . Administrative data — Administrative data can be used to calculate proximity-based parameters and to orient the users geographically. One such global administrative data database is the Global Administrahttp://spatial-analyst.net/worldmaps/ — this repository is constantly updated. http://ngdc.noaa.gov 3 http://srtm.csi.cgiar.org 4 http://topex.ucsd.edu/WWW_html/srtm30_plus.html 5 http://data.gdem.aster.ersdac.or.jp
1 2
99
100
Auxiliary data sources tive Areas (GADM6 ) data set. It comprises borders of countries and lower level subdivisions such as provinces and counties (more than 100,000 areas). Another important global data set is the World Vector Shoreline data set7 at scale 1:250,000 (Soluri and Woodson, 1990). This can be, for example, used to derive the global distance from the sea coast. Eight general purpose thematic layers: boundaries, transportation, drainage, population centers, elevation, vegetation, land use and land cover (al at scale 1:1,000,000) can be obtained via the Global Map Data project8 .
Socio-economic data — The most important global socio-economic data layers are the population density maps and attached socio-economic variables. The Socioeconomic Data and Applications Center (SEDAC9 ) distributes the global population density maps at resolution of 1 km for periods from 1990 up to 2015 (projected density). Water resources — The most detailed and the most accurate inventory of the global water resources is the Global Lakes and Wetlands Database (GLWD10 ), which comprises lakes, reservoirs, rivers, and different wetland types in the form of a global raster map at 30-arcsec resolution (Lehner and Doll, 2004). Shapefiles of the World basins and similar vector data can be best obtained via the Remote Sensing and GIS Unit of the International Water Management Institute (IWMI). Lights at night images — Images of lights at night have shown to be highly correlated with industrial activity and Gross Domestic Product (Doll et al., 2007). A time-series of annual global night light images is available via NOAA’s National Geophysical Data Center11 . The lights at night map contains the lights from cities, towns, and other sites with persistent lighting, including gas flares. The filtered annual composites are available from 1992 until 2003. Land cover maps — Land cover maps are categorical-type maps, commonly derived using semi-automated methods and remotely sensed images as the main input. A Global Land Cover map for the year 2000 (GLC200012 ) at 1 km resolution is distributed by the Joint Research Centre in Italy (Bartholome at al., 2002). A slightly outdated (1998) global map of land cover is the AVHRR Global Land Cover Classification13 , provided at resolutions of 1 and 8 km (Hansen et al., 2000). More detailed land cover maps are distributed nationally. Ellis and Ramankutty (2000) prepared the first global map of the anthropogenic biomes (18 classes) showing dense settlements, villages, croplands, rangelands, forested lands and wildlands. The International Water Management Institute also produced the Global map of Irrigated Areas14 (GMIA; 28 classes) and the Global map of Rainfed Cropped Areas (GMRCA), both at 10 km resolution, and based on twenty years of AVHRR images, augmented with higher resolution SPOT and JERS-1 imagery. Climatic maps — Hijmans et al. (2005) produced global maps of bioclimatic parameters (18) derived (thin plate smoothing splines) using >15,000 weather stations. The climatic parameters include: mean, minimum and maximum temperatures, monthly precipitation and bioclimatic variables15 . Ecoregions / Biogeographic regions — Ecoregions are terrestrial, freshwater and/or marine areas with characteristic combinations of soil and landform that characterize that region. Olson et al. (2001) produced the Terrestrial Ecoregions16 global data set, which shows some 867 distinct eco-units, including the relative richness of terrestrial species by ecoregion. A somewhat more generalized is the FAO’s map of Eco-floristic regions17 (e.g. boreal coniferous forest, tropical rainforest, boreal mountain system etc.). Soil / Geology maps — USGS produced a detailed Global Soil Regions map18 at resolution of 60 arcsec. FAO, IIASA, ISRIC, ISSCAS, JRC have recently produced a 1 km gridded map, merged from various nahttp://biogeo.berkeley.edu/gadm/ http://rimmer.ngdc.noaa.gov/mgg/coast/wvs.html 8 http://www.iscgm.org 9 http://sedac.ciesin.columbia.edu/gpw/ 10 http://www.worldwildlife.org/science/data/item1877.html 11 http://ngdc.noaa.gov/dmsp/global_composites_v2.html 12 http://www-tem.jrc.it/glc2000/ 13 http://glcf.umiacs.umd.edu/data/landcover/data.shtml 14 http://www.iwmigiam.org/info/gmia/ 15 The 1 km resolution maps can be obtained via http://worldclim.org. 16 http://www.worldwildlife.org/science/ecoregions/item1267.html 17 http://cdiac.ornl.gov/ftp/global_carbon/ 18 http://soils.usda.gov/use/worldsoils/mapindex/order.html
6 7
4.1 Global data sets
101
tional soil maps, which is also known as the Harmonized World Soil Database19 (v 1.1). The geological maps are now being integrated via the OneGeology project. USDA Soil Survey Division also distributes a global map of wetlands (includes: upland, lowland, organic, permafrost and salt affected wetlands). International Soil Reference Information Center (ISRIC) maintains a global soil profile database with over 12,000 profiles and over 50 analytical and descriptive parameters (Batjes, 2009). From NOAA’s National Geophysical Data Center one can obtain a point map with all major earth quakes (Significant Earthquake Database20 ; with cca 5000 quakes). Forest / wildlife resources — There are two important global forest/wildlife data sets: (1) The world map of intact forest landscapes21 (hardly touched by mankind) at a scale 1:1,000,000 (includes four classes of intact forests: 1. intact closed forests; 2. intact open forests, 3. woodlands and savannas, closed forests; and 4. open forests, woodlands and savannas) — maintained by the Greenpeace organization (Potapov, P et al., 2008), and (2) World Wilderness Areas22 at scale 1:1,000,000 — distributed via the . UNEP GEO Data Portal (McCloskey and Spalding, 1989). Biodiversity / human impacts maps — Global maps of biodiversity measures for various groups of taxa (e.g. vascular plants, birds and mammals) can be browsed using the World Atlas of Biodiversity viewer (Groombridge and Jenkins, 2002). Similar type of maps can be browsed via the UNEP’s World Conservation Monitoring Centre23 . Kreft and Jetz (2007) recently produced a global map of plant species diversity (number of plant species) by using field records from 1,032 locations. Partners in the GLOBIO consortium created a World Map of Human Impacts24 on the Biosphere. This is basically a map showing a current status of the roads, railways and settlement density. The Carbon Dioxide Information Analysis Center provides access to numerous ecological layers and data of interest to global ecologists. One such product is the Global Biomass Carbon Map25 (Carbon density tones of C ha−1 ), prepared for year 2000 (Ruesch and Gibbs, 2008a).
Fig. 4.1: The 1 km resoluton ETOPO1 Global Relief Model (includes bathymetry data). Elevation is probably the most widely used global environmental data set — it is available globally at resolutions from 1 km up to 90 m (SRTM) i.e. 60 m (GDEM). Elevation can be used to extract over 100 unique land surface parameters.
http://www.fao.org/nr/water/news/soil-db.html http://ngdc.noaa.gov/hazard/earthqk.shtml 21 http://www.intactforests.org/ 22 http://geodata.grid.unep.ch/ 23 http://www.unep-wcmc.org/ 24 http://www.globio.info/region/world/ 25 http://cdiac.ornl.gov/ftp/
19 20
102 4.1.1 Obtaining data via a geo-service
Auxiliary data sources
Much of the geo-data can be accessed directly through i.e. via some of its packages specialized in getting geo-data via some web-service. Paul Wessel, from the University of Hawai’i, maintains A Global Selfconsistent, Hierarchical, High-resolution Shoreline Database26 , which is described in detail in Wessel and Smith (1996). These vectors can be imported to by using the package (see Rgshhs method). Some basic (and slightly out-dated) vector maps are available in the R’s package . This contains data map of political borders, world cities, CIA World Data Bank II27 data, administrative units from the NUTS III (Tertiary Administrative Units of the European Community) and similar. To obtain these vector maps for your GIS, you can run:
R
R
maptools
maps
> > + > + >
library(maps) worldmap <- map2SpatialLines(map("world", fill=TRUE, col="transparent", plot=FALSE), proj4string=CRS("+proj=longlat +ellps=WGS84")) worldmap <- SpatialLinesDataFrame(worldmap, data.frame(name=(map("world"))$names), match.ID=F) writeOGR(worldmap, "worldmap.shp", "worldmap", "ESRI Shapefile")
The original shapefiles of World political borders can be obtained directly via the thematicmapping.org28 . The lower level administrative boundaries (global coverage) can be obtained from the FAO’s GeoNetwork server29 . If you wish to obtain similar type of geographic information but only for a specific point location, you should consider using some of the free web-services such as GeoNames (also available via the package ). For example, to obtain elevation, name of the closest city and/or actual weather at some point location, we can run:
GeoNames
R
> library(geonames) > GNfindNearbyPlaceName(lat=47,lng=9) name lat Atzmännig lng geonameId countryCode countryName fcl fcode distance 47.287633 8.988454 6559633 CH Switzerland P PPL 1.6276
> GNgsrtm3(lat=47,lng=9) 1 srtm3 lng lat 2834 9 47
> GNweather(north=47,east=8,south=46,west=9) clouds weatherCondition 1 few clouds n/a
observation 1 LSZA 231320Z VRB02KT 9999 FEW045 BKN060 04/M02 Q0991 NOSIG ICAO lng temperature dewPoint windSpeed 1 LSZA 8.966667 4 -2 02 humidity stationName datetime lat 1 64 Lugano 2009-01-23 14:20:00 46 hectoPascAltimeter 1 991
Another alternative is Google’s maps service, which allows you to obtain similar information. For example, you can use Google’s geographic services (see also coverage detail of Google maps30 ) to get geographic coordinates given a street + city + country address. First, register your own Google API key. To geocode an address, you can run in :
R
http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html http://www.evl.uic.edu/pape/data/WDB/ 28 http://thematicmapping.org/downloads/world_borders.php 29 http://www.fao.org/geonetwork/srv/en/main.home 30 http://en.wikipedia.org/wiki/Coverage_details_of_Google_Maps
26 27
4.1 Global data sets
103
> readLines(url("http://maps.google.com/maps/geo?q=1600+Amphitheatre+Parkway, + +Mountain+View,+CA&output=csv&key=abcdefg"), n=1, warn=FALSE)
which will give four numbers: 1. HTTP status code, 2. accuracy, 3. latitude, and 4. longitude. In the case from above:
[1] 200.00000 8.00000 37.42197 -122.08414
where the status code is 200 (meaning “No errors occurred; the address was successfully parsed and its geocode has been returned”31 ), the geocoding accuracy is 8 (meaning highly accurate; see also the accuracy constants), longitude is 37.42197 and the latitude is -122.08414. Note that the address of a location needs to be provided in the following format:
"StreetNumber+Street,+City,+Country"
A large number of maps can be also obtained via some of the many commercial WCS’s32 . A popular WMS that allows download of the original vector data is Openstreetmap33 . The original data come in the OSM (Open Street Map) format, but can be easily exported and converted to e.g. ESRI shapefiles using the GIS. Another extensive WMS is NASA’s OnEarth34 .
OpenJUMP
4.1.2
Google Earth/Maps images
You can also consider obtaining color composites of the high resolution imagery (QuickBird, Ikonos) that are 35 36 used in . With package you can automate retrieval and mosaicking of images. For example, to obtain a hybrid satellite image of the Netherlands, it is enough to define the bounding box, position of the center and the zoom level (scale):
Google Earth
RgoogleMaps
> # > + >
library(RgoogleMaps) Get the Maximum zoom level: mzoom <- MaxZoom(latrange=c(50.74995, 53.55488), lonrange=c(3.358871 7.227094), size=c(640, 640))[[1]] mzoom [1] 7
# Get a satellite image of the Netherlands: > MyMap <- GetMap.bbox(center=c(52.1551723, 5.3872035), zoom=mzoom, + destfile="netherlands.png", maptype="hybrid") Read 1 item [1] "http://maps.google.com/staticmap?center=52.15517,5.38720&zoom=7&size=640x640 + &maptype=hybrid&format=png32&key=****&sensor=true" trying URL 'http://maps.google.com/staticmap?center=52.15517,5.38720&zoom=7 + &size=640x640&maptype=hybrid&format=png32&key=****=true' Content type 'image/png' length 703541 bytes (687 Kb) opened URL downloaded 687 Kb netherlands.png has GDAL driver PNG and has 640 rows and 640 columns > PlotOnStaticMap(MyMap, lat=52.1551723, lon=5.3872035)
See also the status code table at: http://code.google.com/apis/maps/documentation/reference.html. http://www.ogcnetwork.net/servicelist — a list of Open Geospatial Consortium (OGC) WMS’s. 33 http://www.openstreetmap.org/ — see export tab. 34 http://onearth.jpl.nasa.gov/ 35 This has some copyright restrictions; see http://www.google.com/permissions/geoguidelines.html. 36 http://cran.r-project.org/web/packages/RgoogleMaps/
31 32
104
Auxiliary data sources
The tiles obtained using are blocks of maximum 640×640 pixels distributed as PNG or JPG (compressed) images, which makes them of limited use compare to the multi-band satellite images we receive from the original distributors. Nevertheless, Google contains high resolution imagery of high spatial quality (Potere, 2008), which can be used to extract additional content for a smaller size GIS e.g. to digitize forest borders, stream lines and water bodies. Before you can extract content, you need to attach coordinates to the Static map i.e. you need to georeference it. To achieve this, we first need to estimate the bounding box coordinates of the image. This is possible with the help of the XY2LatLon.R script:
RgoogleMaps
# > + > > # > + > >
Get the XY2LatLon.R script from CRAN: download.file("http://cran.r-project.org/src/contrib/RgoogleMaps_1.1.6.tar.gz", destfile=paste(getwd(), "/", "RgoogleMaps.tar.gz", sep="")) library(R.utils) gunzip("RgoogleMaps.tar.gz", overwrite=TRUE, remove=TRUE) under windows, you need to use 7z software to unzip *.tar archives: download.file("http://downloads.sourceforge.net/sevenzip/7za465.zip", destfile=paste(getwd(), "/", "7za465.zip", sep="")) unzip("7za465.zip") system("7za e -ttar RgoogleMaps.tar XY2LatLon.R -r -aos") 7-Zip (A) 4.65 Copyright (c) 1999-2009 Igor Pavlov Processing archive: RgoogleMaps.tar Extracting RgoogleMaps\R\XY2LatLon.R 2009-02-03
Everything is Ok Size: 1001 Compressed: 1228800 > source("XY2LatLon.R") # Read the Google Static image into R: > png.map <- readGDAL("netherlands.png") netherlands.png has GDAL driver PNG and has 640 rows and 640 columns # > + > + > > # > > Estimate the bounding box coordinates: bbox.MyMap <- data.frame(XY2LatLon(MyMap, X=png.map@bbox[1,]-640/2, Y=png.map@bbox[2,]-640/2)) google.prj <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" coordinates(bbox.MyMap) <- ∼ lon+lat proj4string(bbox.MyMap) <- CRS("+proj=longlat +ellps=WGS84") Estimate coordinates in Google Maps projection system: bbox.google.prj <- spTransform(bbox.MyMap, CRS(google.prj)) bbox.google.prj SpatialPoints: lon lat [1,] 756387.0 7035766 [2,] 768616.9 7047996 Coordinate Reference System (CRS) arguments: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
then attach the right georeference (east, west, north, south coordinates):
# copy the bounding box coordinates: > png.map@bbox[1,1] <- bbox.google.prj@coords[1,1]
4.1 Global data sets
105
> > > # > + > > > > # > + # > # > > +
png.map@bbox[2,1] <- bbox.google.prj@coords[1,2] png.map@bbox[1,2] <- bbox.google.prj@coords[2,1] png.map@bbox[2,2] <- bbox.google.prj@coords[2,2] cell size: png.map@grid@cellsize <- round(c((png.map@bbox[1,2]-png.map@bbox[1,1])/640, (png.map@bbox[2,2]-png.map@bbox[2,1])/640), 1) png.map@coords[1,1] <- png.map@bbox[1,1]+png.map@grid@cellsize[1]/2 png.map@coords[1,2] <- png.map@bbox[2,1]+png.map@grid@cellsize[2]/2 png.map@coords[2,1] <- png.map@bbox[1,2]-png.map@grid@cellsize[1]/2 png.map@coords[2,2] <- png.map@bbox[2,2]-png.map@grid@cellsize[2]/2 cell offset: png.map@[email protected] <- c(png.map@bbox[1,1]+png.map@grid@cellsize[1]/2, png.map@bbox[2,1]+png.map@grid@cellsize[2]/2) attach the correct prj: proj4string(png.map) <- CRS(google.prj) Export map to GIS format: write.asciigrid(png.map[1], "netherlands_B1.asc", na.value=-1) writeGDAL(netherlands[c(1,2,3)], "netherlands.tif", drivername="GTiff", type="Byte", options="INTERLEAVE=PIXEL")
Note that Google Maps currently uses a modification of the Mercator projection system37 to display tiles. To reproject this grid to some other system consider using the proj4 functionality, as explained further in §5.6.2. We can also just quickly estimate how many tiles were used to generate this map by using the maptiler script38 (courtesy of Klokan Petr Pˆidal): r
SAGA
Python
# You need to first install and register Python on your machine! > download.file("http://www.maptiler.org/google-maps-coordinates-tile-... [TRUNCATED] + destfile=paste(getwd(), "/", "globalmaptil ..." ... [TRUNCATED] trying URL 'http://www.maptiler.org/google-...-projection/globalmaptiles.py' Content type 'text/plain' length 16529 bytes (16 Kb) opened URL downloaded 16 Kb > system(paste("python.exe globalmaptiles.py", mzoom, NLprovs.ll@bbox[2,1], + NLprovs.ll@bbox[1,1], NLprovs.ll@bbox[2,2], NLprovs.ll@bbox[1,2])) Spherical Mercator (ESPG:900913) coordinates for lat/lon: (373907.80392051308, 6577181.3183709756) Spherical Mercator (ESPG:900913) cooridnate for maxlat/maxlon: (804516.38277077128, 7086303.4059718344) 7/65/85 ( TileMapService: z / x / y ) Google: 65 42 Quadkey: 1202021 ( 6281 ) EPSG:900913 Extent: (313086.06785608083, 6574807.4249777198, 626172.13571216539, 6887893.4928338043) WGS84 Extent: (50.736455137010637, 2.8124999999999902, 52.482780222078226, 5.6250000000000133) gdalwarp -ts 256 256 -te 313086.067856 6574807.42498 626172.135712 6887893.49283 <your-raster-file-in-epsg900913.ext> 7_65_85.tif 7/66/85 ( TileMapService: z / x / y ) Google: 66 42 Quadkey: 1202030 ( 6284 ) EPSG:900913 Extent: (626172.13571216539, 6574807.4249777198,
37 38
http://www.spatialreference.org/ref/user/google-projection/ http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/
106
Auxiliary data sources
939258.20356824622, 6887893.4928338043) WGS84 Extent: (50.736455137010637, 5.6250000000000133, 52.482780222078226, 8.4375000000000036) gdalwarp -ts 256 256 -te 626172.135712 6574807.42498 939258.203568 6887893.49283 <your-raster-file-in-epsg900913.ext> 7_66_85.tif 7/65/86 ( TileMapService: z / x / y ) Google: 65 41 Quadkey: 1202003 ( 6275 ) EPSG:900913 Extent: (313086.06785608083, 6887893.4928338043, 626172.13571216539, 7200979.5606898852) WGS84 Extent: (52.482780222078226, 2.8124999999999902, 54.162433968067809, 5.6250000000000133) gdalwarp -ts 256 256 -te 313086.067856 6887893.49283 626172.135712 7200979.56069 <your-raster-file-in-epsg900913.ext> 7_65_86.tif 7/66/86 ( TileMapService: z / x / y ) Google: 66 41 Quadkey: 1202012 ( 6278 ) EPSG:900913 Extent: (626172.13571216539, 6887893.4928338043, 939258.20356824622, 7200979.5606898852) WGS84 Extent: (52.482780222078226, 5.6250000000000133, 54.162433968067809, 8.4375000000000036) gdalwarp -ts 256 256 -te 626172.135712 6887893.49283 939258.203568 7200979.56069 <your-raster-file-in-epsg900913.ext> 7_66_86.tif
which shows that four tiles are needed to represent the whole of the Netherlands at zoom level 7. Again, you need to be aware that these maps/images are copyrighted, so you should really use them only for your personal purpose. In addition, access to Google imagery is possible only via the Google maps API, which means that you first need to register your key etc. If you really require detailed remotely sensed images for larger study area, then consider obtaining some of the popular RS products (Landsat, Spot, ASTER, Ikonos) from the official data distributers. A number of similar Python scripts can be found at GDAL utilities section39 . For example, if you have a relatively large raster map and you want to export it to KML (Super Overlays), you can create a type of hierarchical images by using the gdal2tiles.py script (also available through ). This will automatically split the map into tiles, generate a directory with small tiles and metadata following the OSGeo Tile Map Service Specification40 .
Google Earth FWtools
4.1.3 Remotely sensed images Remotely sensed images are increasingly the main source of data for many national and continental scale mapping projects. The amount of field and remotely sensed data in the world is rapidly increasing. To get an idea about how many sensors are currently available, see for example the ITC’s database41 . The most common satellites/images with global coverage that are available at low cost or for free, and which are of interest for global modeling projects, are (Fig. 4.2): Landsat Landsat collects multispectral satellite imagery at resolutions 15–30 meters. A number of their products is also available at no cost. High resolution (15 m) Landsat images for nearly all of the world (for years 1990 and 2000) can be downloaded from the NASA’s Seamless Server42 . European (harmonized) mosaic of Landsat images is distributed by JRC Ispra (see Image200043 ). Another excellent repository of free global imagery is the GLCF geo-portal44 operated by the University of Maryland.
http://www.gdal.org/gdal_utilities.html http://code.google.com/apis/kml/articles/raster.html 41 http://www.itc.nl/research/products/sensordb/AllSensors.aspx 42 http://seamless.usgs.gov/ 43 http://image2000.jrc.ec.europa.eu 44 http://glcf.umiacs.umd.edu/portal/geocover/
39 40
4.1 Global data sets
107
SPOT SPOT is a commercial distributor of high quality resolution satellite imagery (multispectral images at resolution of 2.5–20 m). The most recent satellite SPOT 6 carries a High Resolution Stereoscopy (HRS) sensor that allows production of 3D imagery. SPOT vegetation45 offers relatively coarse vegetationbased 10–day images of the whole Earth collected in the period from 1998 until today. Only two bands are available at the moment: NDVI and radiometry images. Ikonos The IKONOS sensor (Satellite) is a high-resolution commercial satellite operated by GeoEye company. The images are available at two resolutions: 3.2 m (multispectral, Near-Infrared), and 0.82 m (panchromatic) (Dial et al., 2003). Ikonos images are sold per km2 (a standard scene is of size 12,5×12,5 km); the user defines an area of interest of any shape but with a minimum surface area. Archived Ikonos images at discounted price can be obtained via various companies. Meteosat The Meteosat46 Second Generation (MSG) satellites (from Meteosat-8 onwards) produce SEVIRI (Spinning Enhanced Visible and Infrared Imagery) 15–minutes meteorological images at a resolution of 1 km. The most attractive data set for environmental applications is the High Rate SEVIRI, which consists of 12 spectral channels including: visible and near infrared light, water vapor band, carbon dioxide and ozone bands. ENVISAT The ENVISAT satellite is a platform for several instruments adjusted for monitoring of the environmental resources: ASAR, MERIS, AATSR, MWR and similar. The MEdium Resolution Image Spectrometer (MERIS47 ) is used to obtain images of the Earth’s surface at a temporal resolution of 3–days. The images comprise 15 bands, all at a resolution of 300 m. To obtain MERIS images (Category 1 use data), one needs to register and wait few days to receive access to the repository48 (unlike the MODIS images that are available directly via an FTP). AVHRR Although outdated, AVHRR was the one of the main global environmental monitoring systems in the 80’s. Principal components derived from a set of 232 monthly NDVI images49 can be obtained from the author’s repository of world maps. The original images are available at a resolution of 300 arcseconds (cca 10 km), and cover the period July, 1981 through September, 2001. MODIS MODIS contains a number of products ranging from raw multispectral images, to various vegetation and atmospheric indices at a resolution of 250 m (also available at coarser resolutions of 500 m and 1 km), and at very high temporal resolution. If you only wish to use completed MODIS products, then look at NASA’s Earth Observation (NEO50 ) portal, which distributes global time-series of MODISderived parameters such as: snow cover and ice extent, leaf area index, land cover classes, net primary production, chlorophyll concentration in the sea and sea surface temperature, cloud water content, carbon monoxide in the atmosphere and many more. All global maps on NEO are freely available for public use. You can simply download the 0.1 arcdegrees (∼10 km) GeoTiff’s and load them to your GIS. Please credit NASA as the source. From the RS systems listed above, one needs to be emphasized — NASA’s MODIS Earth observation system — possibly one of the richest sources of remote-sensing data for monitoring of environmental dynamics (Neteler, 2005; Lunetta et al., 2006; Ozdogana and Gutman, 2008). This is due to the following reasons: (1.) it has global coverage; (2.) it has a relatively high temporal resolution/coverage (1–2 days; Fig. 4.2); (3.) it is open-access (see the MODIS licence specification); (4.) a significant work has been done to filter the original raw images for clouds and artifacts; a variety of complete MODIS products such as composite 15–day and monthly images is available at three resolutions: 250 m, 500 m and 1 km;
http://www.spot-vegetation.com http://www.eumetsat.int 47 http://envisat.esa.int/instruments/ 48 http://earth.esa.int/dataproducts/accessingeodata/ 49 Available via the International Water Management Institute at http://iwmidsp.org/. 50 http://neo.sci.gsfc.nasa.gov/
45 46
108
Auxiliary data sources
Global environmental monitoring
n ch Te og ol y 08 20
gy lo no ch Te
gy lo no ch Te
00 20
15 20
LANDSAT-7 Revisit (days)
SPOT-4, 5 IKONOS-2 IRS-P6
ENVISAT SPOT-6
Urban planning Civil engineering Land use mapping Forestry, coastal resources
MODIS GEOSAT
Fig. 4.2: Resolution and revisit time of some common imaging satellites. Modified after Davis et al. (2009).
(5.) efficient tools exist to obtain various MODIS products and import them to various GIS; One of the best known MODIS products51 for terrestrial environmental applications is the Enhanced Vegetation Index (EVI), which is the improved NDVI (Huete et al., 2002). EVI corrects distortions in the reflected light caused by particles in the air as well as ground cover below the vegetation. The EVI also does not become saturated as easily as the NDVI when viewing rainforests and other areas with large amounts of chlorophyll. EVI can be directly related to the photosynthetic production of plants, and indirectly to the green biomass (Huete et al., 2002). By observing dynamics of EVI for an area of the Earth’s surface, we can conclude about the vegetation dynamics within a season, but also detect long-term trends and sudden changes in the biomass (e.g. due to forest fires, deforestation, urban growth and similar).
4.2
Download and preparation of MODIS images
This section explains how to automate download, mosaicking, resampling and import of MODIS product into a GIS. We focus on the Land Surface Temperature (LST) images, which are subsequently used to improve spatio-temporal interpolation of temperatures in chapter 11. Before you can start downloading MODIS images from within , you need to obtain and install some necessary applications and libraries:
R
R
RCurl52 — allows you to list directories on an FTP server; wget53 — this will automate download of images from within R; simply put the wget.exe in your
Windows system folder54 ;
MRT55 — MODIS reproject tool can be used to mosaic MODIS images, resample them to other coordinate systems, and export images to more common GIS formats; After you have finished installing all these software programs, you also need to specify the location of MRT and the directory where you want to output all EVI maps you will produce:
https://lpdaac.usgs.gov/lpdaac/products/modis_product_table http://www.omegahat.org/RCurl/ 53 http://users.ugent.be/~bpuype/wget/
51 52 55
54 Note: make sure you disable your antivirus tools such as Norton or McAfee otherwise it might block has not yet been tested under Mac OS X.
wget from running. This script
https://lpdaac.usgs.gov/lpdaac/tools/modis_reprojection_tool
4.2 Download and preparation of MODIS images
109
> > # > # > # > >
library(RCurl) library(rgdal) location of the mosiacing tool: MRT <- 'E:\\MODIS\\MRT\\bin\\' location of the working directory: workd <- 'E:\\MODIS\\HR\\' location of the MODIS 1 km blocks: MOD11A2 <- "ftp://e4ftl01u.ecs.nasa.gov/MOLT/MOD11A2.005/" MOD11A2a <- "ftp://anonymous:[email protected]/MOLT/MOD11A2.005/"
MODIS images are typically distributed as HDF (Hierarchical Data Format) 10 by 10 arcdegree-tiles, projected in the sinusoidal projection (Fig. 4.3). The sinusoidal projection has been promoted by geographers as the most suited projection for global image databases (Chang Seong et al., 2002). Unfortunately, both HDF format and sinusoidal projection are yet not supported in many GIS. Therefore, before you can use MODIS images, you will need to run some pre-processing to glue the tiles and convert the data into a more usable format. We wish to obtain the 8–day 1 km resolution MODIS LST images56 , i.e. the average values of clear-sky LSTs during the 8–day periods. We will first download the original tiles, then mosaic them and resample them to the UTM projection system (zone 33), and then export them to a more common GIS format (GeoTiff). The MOD13A3 tiles can be browsed and downloaded directly via NASA’s FTP server. However, for Fig. 4.3: MODIS HDF tiles (Sinusoidal grid). larger areas, you will often need to download tens of tiles, so that is fairly useful to automate the processing. We start by fetching the list of sub-directories (the dates) of interest:
# get the list of directories: > items <- strsplit(getURL(MOD11A2), "\n")[[1]] > items[2] [1] "drwxr-xr-x 2 90 129024 Dec 12 2008 2000.03.05\r" # # > # > > > you get the folder names but in form of a unix directory listing get the last word of any lines that start with 'd': folderLines <- items[substr(items, 1, 1)=='d'] get the directory names and create a new data frame: dirs <- unlist(lapply(strsplit(folderLines, " "), function(x){x[length(x)]})) dates <- data.frame(dirname=unlist(strsplit(dirs, "\r"))) str(dates)
'data.frame':
430 obs. of 1 variable: $ dirname: Factor w/ 430 levels "2000.03.05","2000.03.13",..: 1 2 3 4 5 6 ...
which gives 430 dates; we are only interested in year 2006 (i.e. 268–313; i.e. 45 dates). Note that there are some differences between MS-Windows vs Mac OS machines in the use of "\r\n" vs "\r" ("\n" is echoed out in the terminal with the file names). Next, we need to known the h/v position of the MODIS blocks. For example, we want to generate EVI images for the whole area of Croatia. This area covers two MODIS tiles: h18v04 and h19v04 (Fig. 4.3). Each MODIS tile has a unique name. We can do a directory listing and get the full names of the tiles on the FTP by combining the getURL and grep methods:
R
56
This data set is known by the name MOD11A2.
110
Auxiliary data sources
> getlist <- strsplit(getURL(paste(MOD11A2, dates$dirname[[1]], "/", sep=""), + .opts=curlOptions(ftplistonly=TRUE)), "\r\n")[[1]] > str(getlist) chr [1:1268] "BROWSE.MOD11A2.A2006001.h00v08.005.2008098013929.1.jpg" ...
This means that the FTP directory of interest57 contains a total of 1268 files. We wish to obtain only the names of the HDF files (two tiles) for our area of interest. These can be obtained using the grep method:
> > > + > + >
dates$BLOCK1 <- rep(NA, length(dates$dirname)) dates$BLOCK2 <- rep(NA, length(dates$dirname)) BLOCK1 <- getlist[grep(getlist, pattern="MOD11A2.********.h18v04.*************.hdf")[1]] BLOCK2 <- getlist[grep(getlist, pattern="MOD11A2.********.h19v04.*************.hdf")[1]] BLOCK1 [1] "MOD11A2.A2006001.h18v04.005.2008098031505.hdf"
which means that we have successfully determined names of the tiles we are interested in downloading. We look only for the first (HDF) file; the second file with the same name but .XML extension carries the production metadata58 . Next, we can download each tile using the download.file method, with help from the package59 :
wget
> download.file(paste(MOD11A2a, dates$dirname[[1]], "/", BLOCK1,sep=""), + destfile=paste(getwd(), "/", BLOCK1, sep=""), mode='wb', method='wget') --2009-03-01 18:21:37-- ftp://anonymous:*password*@e4ftl01u.ecs.nasa.gov/MOLT/... => `D:/MODIS/HR/MOD11A2.A2006001.h18v04.005.2008098031505.hdf' Resolving e4ftl01u.ecs.nasa.gov... 152.61.4.83 Connecting to e4ftl01u.ecs.nasa.gov|152.61.4.83|:21... connected. Logging in as anonymous ... Logged in! ==> SYST ... done. ==> PWD ... done. ==> TYPE I ... done. ==> CWD /MOLT/MOD11A2.005/2006.01.01 ... done. ==> SIZE MOD11A2.A2006001.h18v04.005.2008098031505.hdf ... 6933421 ==> PASV ... done. ==> RETR MOD11A2.A2006001.h18v04.005.2008098031505.hdf ... done. Length: 6933421 (6.61M)
Once we have downloaded all tiles of interest, we can mosaic them into a single image. To do this, we use the MODIS Resampling Tool, that you should have already installed on your machine. We first generate a parameter file containing a list of tiles that need to be mosaicked:
> > > >
mosaicname <- file(paste(MRT, "TmpMosaic.prm", sep=""), open="wt") write(paste(workd, BLOCK1, sep=""), mosaicname) write(paste(workd, BLOCK2, sep=""), mosaicname, append=T) close(mosaicname)
and then run the MRT mosaicking tool using this parameter file:
# Generate a mosaic: > shell(cmd=paste(MRT, 'mrtmosaic -i ', MRT, 'TmpMosaic.prm -s + "1 0 0 0 0 0 0 0 0 0 0 0" -o ', workd, 'TmpMosaic.hdf', sep="")) ****************************************************************************** MODIS Mosaic Tool (v4.0 February 2008) Start Time: Thu Jul 30 14:05:26 2009
57 58 59
ftp://e4ftl01u.ecs.nasa.gov/MOLT/MOD11A2.005/2006.01.01/
Much about the HDF tile can be read already from the file name; see also the MODIS Naming Conventions. You need to download the wget.exe to your windows system directory, otherwise you will not be able to download the tiles from R.
4.2 Download and preparation of MODIS images
111
-----------------------------------------------------------------Input filenames (2): D:\MODIS\HR\MOD11A2.A2006001.h18v04.005.2008098031505.hdf D:\MODIS\HR\MOD11A2.A2006001.h19v04.005.2008098031444.hdf Output filename: D:\MODIS\HR\TmpMosaic.hdf Mosaic Array: file[ 0] file[ 1] Mosaic : processing band LST_Day_1km % complete (1200 rows): 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% Output mosaic image info -----------------------output image corners (lat/lon): UL: 50.000000000003 0.000000000000 UR: 50.000000000003 31.114476537210 LL: 39.999999999999 0.000000000000 LR: 39.999999999999 26.108145786645 output image corners (X-Y projection units): UL: 0.000000000000 5559752.598833000287 UR: 2223901.039532999974 5559752.598833000287 LL: 0.000000000000 4447802.079065999947 LR: 2223901.039532999974 4447802.079065999947 band 1) LST_Day_1km type lines smpls pixsiz UINT16 1200 2400 926.6254 min 7500 max 65535 fill 0
End Time: Thu Jul 30 14:05:28 2009 Finished mosaicking! ******************************************************************************
where "-s 1 0 0 0 0 0 0 0 0 0 0" is the definition of the spectral selection (1st band is the 8–Day daytime 1 km grid land surface temperature), and 'TmpMosaic.hdf' is the temporary mosaic HDF image. Next, we need to generate a MRT parameter file that can be use to resample the image to the target coordinate system:
# > > > > > > > > + > + > > + > > > > > > > >
resample to UTM 33N: filename <- file(paste(MRT, "mrt2006_01_01.prm", sep=""), open="wt") write(paste('INPUT_FILENAME = ', workd, 'TmpMosaic.hdf', sep=""), filename) write(' ', filename, append=TRUE) write('SPECTRAL_SUBSET = ( 1 )', filename, append=TRUE) write(' ', filename, append=TRUE) write('SPATIAL_SUBSET_TYPE = OUTPUT_PROJ_COORDS', filename, append=TRUE) write(' ', filename, append=TRUE) write(paste('SPATIAL_SUBSET_UL_CORNER = (', XUL, YUL, ')'), filename, append=TRUE) write(paste('SPATIAL_SUBSET_LR_CORNER = (', XLR, YLR, ')'), filename, append=TRUE) write(' ', filename, append=TRUE) write(paste('OUTPUT_FILENAME = ', workd, 'LST', dirname1, '.tif', sep=""), filename, append=TRUE) write(' ', filename, append=TRUE) write('RESAMPLING_TYPE = NEAREST_NEIGHBOR', filename, append=TRUE) write(' ', filename, append=TRUE) write('OUTPUT_PROJECTION_TYPE = UTM', filename, append=TRUE) write(' ', filename, append=TRUE) write('OUTPUT_PROJECTION_PARAMETERS = ( ', filename, append=TRUE) write(paste(lon.c, lat.c, '0.0'), filename, append=TRUE) write(' 0 0.0 0.0', filename, append=TRUE)
112
Auxiliary data sources
> > > > > > > > > >
write(' 0.0 0.0 0.0', filename, append=TRUE) write(' 0.0 0.0 0.0', filename, append=TRUE) write(' 0.0 0.0 0.0 )', filename, append=TRUE) write(' ', filename, append=TRUE) write('UTM_ZONE = 33', filename, append=TRUE) write('DATUM = WGS84', filename, append=TRUE) write(' ', filename, append=TRUE) write('OUTPUT_PIXEL_SIZE = 1000', filename, append=TRUE) write(' ', filename, append=TRUE) close(filename)
and we can again run the MRT to resample the map:
> shell(cmd=paste(MRT, 'resample -p ', MRT, 'mrt2006_01_01.prm', sep="")) ****************************************************************************** MODIS Reprojection Tool (v4.0 February 2008) Start Time: Thu Jul 30 14:09:29 2009 -----------------------------------------------------------------Input image and reprojection info --------------------------------input_filename: D:\MODIS\HR\TmpMosaic.hdf output_filename: D:\MODIS\HR\LST2006_01_01.tif input_filetype: HDF-EOS output_filetype: GEOTIFF input_projection_type: SIN input_datum: WGS84 output_projection_type: UTM output_zone_code: 33 output_datum: WGS84 resampling_type: NN input projection parameters: 6371007.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 86400.00 0.00 0.00 0.00 0.00 0.00 0.00 output projection parameters: 16.33 44.35 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 input image corners (lat/lon): UL: 50.00 0.00 UR: 50.00 31.11 LL: 40.00 0.00 LR: 40.00 26.11 input image spatial subset corners (lat/lon): UL: 46.55 13.18 UR: 46.47 19.55 LL: 42.17 13.31 LR: 42.11 19.23 band 1) LST_Day_1km select type lines smpls pixsiz 1 UINT16 1200 2400 926.6254 min max 7500 65535 fill 0
SINUSOIDAL PROJECTION PARAMETERS: Radius of Sphere: 6371007.181000 meters Longitude of Center: 0.000000 degrees False Easting: 0.000000 meters False Northing: 0.000000 meters
4.2 Download and preparation of MODIS images
113
UNIVERSAL TRANSVERSE MERCATOR (UTM) PROJECTION PARAMETERS: Zone: 33 Semi-Major Axis of Ellipsoid: Semi-Minor Axis of Ellipsoid: Scale Factor at C. Meridian: Longitude of Central Meridian: 6378137.000000 meters 6356752.314245 meters 0.999600 15.000000 degrees
NNResample : processing band LST_Day_1km % complete (487 rows): 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% Output image info ----------------output image extents (lat/lon): UL: 46.548860613555 13.175657895706 UR: 46.472377855130 19.561116709234 LL: 42.166491811106 13.306850810151 LR: 42.100790881433 19.234156734254 output image extents (X-Y projection units): UL: 360141.478220875026 5156649.305139659904 UR: 850141.478220875026 5156649.305139659904 LL: 360141.478220875026 4669649.305139659904 LR: 850141.478220875026 4669649.305139659904 band 1) LST_Day_1km type lines smpls pixsiz UINT16 487 490 1000.0000 min 7500 max fill 65535 0
End Time: Thu Jul 30 14:09:30 2009 Finished processing! ******************************************************************************
Fig. 4.4: A sample of downloaded and resampled MODIS LST images showing the average values of clear-sky land surface temperature (C°) during an 8–day period; see further chapter 11.
By putting these operations in a loop, one can automate downloading, mosaicking and resampling of MODIS images for large areas. An example of such a script can be obtained from the book’s homepage. Just to check that everything is OK with the maps, we can use the GDALinfo (a plot of a sample of images is shown in Fig. 4.4):
> GDALinfo("LST2006_01_01.LST_Day_1km.tif") rows columns 487 490
114
Auxiliary data sources
bands origin.x origin.y res.x res.y oblique.x oblique.y driver projection file
1 360641.5 4669149 1000 1000 0 0 GTiff +proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs LST2006_01_01.LST_Day_1km.tif
4.3 Summary points
The availability of remotely sensed data is everyday increasing. Every few years new generation satellites are launched that produce images at finer and finer detail, at shorter and shorter revisit time, and with richer and richer content (Fig. 4.2). For example, DEMs are now available from a number of sources. Detailed and accurate images of topography can now be ordered from remote sensing systems such as SPOT and ASTER; SPOT5 offers the High Resolution Stereoscopic (HRS) scanner, which can be used to produce DEMs at resolutions of up to 5 m. The cost of data is either free or dropping in price as technology advances. Likewise, access to remotely sensed data and various thematic maps is becoming less and less a technical problem — live image archives and ordering tools are now available. Multi-source layers are also increasingly compatible considering the coverage, scale/resolution, format and metadata description. The intergovernmental Group of Earth Observation is building the cyber-infrastructure called GEOSS60 needed to enhance merging of various geoinformation layers to a single multi-thematic, multi-purpose data repository. From the time when satellite sensors were used for military purposes only, we now live in an era when anybody can load high resolution images of earth and detailed digital maps to their computer. In the context of geostatistical mapping, the number of possible covariates is today increasing dramatically, so that there is practically no reason any more to study point-sampled variables in relation to its spatial coordinates only (Pebesma, 2006). External variables should be used wherever possible to improve geostatistical mapping. Majority of global spatial layers are today available at no cost at resolutions 1–10 km (see the world maps repository); MODIS imagery is freely available at resolution of 250 m; global elevation data is available at resolutions of 60–90 m; archive Landsat images (15–30 m) can be obtained for most of the world. All this proves that there is no reason not to use auxiliary covariates (i.e. methods such as regression-kriging) for geostatistical mapping. In this chapter, I have specifically put an emphasis on MODIS products, mainly because the images are relatively easy to obtain for various parts of the world (or globally), and because MODIS images are valuable for monitoring ecosystem dynamics. From the aspect of ecology and nature conservation, the importance of MODIS images for environmental management and nature conservation is enormous (Lunetta et al., 2006). They can be obtained by anybody at any time and can be used as an independent source of information to quantify degradation of natural systems that is possibly not visible using local sources of information. MODIS EVI images can be used to show changes in land cover, deforestation, damage caused by global warming or even fine long-term succession processes. Because MODIS EVI is available from year the 2000, anybody can compare the current situation with the situation from a decade ago. The script presented can be obtained form the book’s homepage and adopted to download any type of MODIS products available from the USGS Land Processes Distributed Active Archive Center. The next generation remote sensing systems/surveys will certainly be multi-thematic and based on active surface-penetrating sensors. The future of automated mapping lays in using technologies such as LiDAR in combination with other optical, radar and hyperspectral sensors. This type of multi-thematic imagery will enable analysts to represent both surface and sub-surface properties of objects of interest, so that none of the important characteristics are overlooked. Multi-thematic scanners should also make the airborne surveys cheaper and hence suited for local and regional scale environmental studies.
60
http://www.earthobservations.org
4.3 Summary points
115
Further reading: Doll, C.N.H., Muller, J.-P Morley, J.G. 2007. Mapping regional economic activity from night-time light ., satellite imagery. Ecological Economics, 57(1): 75–92. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P .G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965–1978. Pebesma, E.J., 2006. The Role of External Variables and GIS Databases in Geostatistical Analysis. Transactions in GIS, 10(4): 615–632.
http://www.fao.org/geonetwork/srv/en/main.home — FAO’s GeoNetwork server. http://www.geoportal.org/web/guest/geo_image_gallery — ESA’s GeoPortal gallery. http://geodata.grid.unep.ch/ — UNEP/GRID GEO DataPortal. http://neo.sci.gsfc.nasa.gov/ — NASA’s Earth Observation (NEO) portal. http://glcf.umiacs.umd.edu/portal/geocover/ — Global Land Cover Facility (GLCF) geoportal operated by the University of Maryland.
https://lpdaac.usgs.gov/lpdaac/ — Land Processes Distributed Active Archive Center (distributor of MODIS products).