I am tring to read the following GIS datafile into R:
Some Data Specifications are here:
Image Type: Generic Flat Binary, Byte Interleave By Line (BIL)
Projection: Interrupted Goode Homolosine
I tried to use the R package "raster", but failed.
library(raster)
r <- raster(file.choose())
Error in .local(.Object, ...) :
`C:\global_forest_cover.img' not recognised as a supported file format.
Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", :
Cannot create a RasterLayer object from this file.
So, how should I load this GIS data into R? Also, convert the projection "Interrupted Goode Homolosine" to LongLat ?
Create a file called global_forest_cover.hdr
and stick the header info from http://edc2.usgs.gov/glcc/fao/header_file.php into it:
BYTEORDER M
LAYOUT BIL
NROWS 15059
NCOLS 36543
NBANDS 1
NBITS 8
BANDROWBYTES 36543
TOTALROWBYTES 36543
ULXMAP -17359000
ULYMAP 8673000
XDIM 1000
YDIM 1000
then read the .img
file:
> forest = raster("global_forest_cover.img")
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.6.3, released 2009/11/19
Path to GDAL shared files: /usr/share/gdal16
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files: (autodetected)
> forest
class : RasterLayer
dimensions : 15059, 36543, 550301037 (nrow, ncol, ncell)
resolution : 1000, 1000 (x, y)
extent : -17359500, 19183500, -6385500, 8673500 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/nobackup/rowlings/Downloads/global_forest_cover.img
names : global_forest_cover
values : 0, 255 (min, max)
GDAL sees the .hdr
file and uses that to work out the structure of the 15059*36543 bytes in the .img
file.
However it does not have the projection information. But even if it did, warping an IGH projection back to lat-long might be somewhat problematic. If you look at other IGH projected maps you'll see they involve tearing down the oceans to help flatten the globe. To get back to lat-long you need to reverse all these flattenings and tears. It should be possible, and I think the latest PROJ4 libraries have support for IGH projections, but maybe only in one direction. BUT unless you use exactly the same code that was used to create the data you have but in reverse you might not get the right answer out.
Given that the .hdr
file wasn't included in the zip, and that the projection isn't really a good one for further analysis, I'd go looking for another source of the data. This data is clearly only meant for display purposes. And its also poor because it doesn't seem to have missing data values where the IGH projection has split the earth - compare the pictures on the http://edc2.usgs.gov/glcc/fao/index.php page with other IGH projections http://en.wikipedia.org/wiki/Goode_homolosine_projection
I thought the data here might be more appropriate: http://www.fao.org/forestry/32203/en/
> r=raster("./fceurope/europe/w001001.adf")
> r
class : RasterLayer
dimensions : 20000, 40000, 8e+08 (nrow, ncol, ncell)
resolution : 0.009, 0.009 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
data source : /home/nobackup/rowlings/Downloads/fceurope/europe/w001001.adf
names : w001001
values : 1, 6 (min, max)
Raster Attribute Table
fields : ID COUNT
min : 1 129626
max : 6 4963775
But I'm slightly confused between the 'europe' in the name and the global extent. Its also a very very large raster. Ah, its got Europe in the middle of nowhere. Note the projection is now lat-long. You may want to thin the grids out before putting them together to create a global raster from all the files given there.