[GIS] How to read *.adf files into R

dataesri-grid-formatrraster

I want to load .adf files into R. The data is from this page:
http://www.fao.org/geonetwork/srv/en/metadata.show?id=14057

I tried the following code that I found after some research in the internet.
The problem is, that in the class RasterLayer I get negativ values that shouldn't be there. I don't know why this happens, so hopefully someone can help me!?

Code:

library(rgdal)
library(RColorBrewer)
dpath<- path...

x <- new("GDALReadOnlyDataset", dpath)
getDriver(x)
getDriverLongName(getDriver(x))
xx<-asSGDF_GROD(x)
r <- raster(xx)

The output for 'r' is:

r
class : RasterLayer
dimensions : 2160, 4320, 9331200 (nrow, ncol, ncell)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : in memory
names : band1
values : -997, 16 (min, max)

The '16' in the values refers to the 16 classes of length of growing period. But I wonder where those '-997' come from. Maybe something wrong with the coord. ref?

Here is also a data summary of 'xx':

Data summary:
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-997 3 5 -9 8 16 7123158

And if we look at the data in xx more closely:

table(xx$band1)

-997 1 2 3 4 5 6 7 8 9 10 11 12
31711 429643 83011 166674 207228 270161 240958 183342 118608 98795 88473 73743 56022
13 14 15 16
30104 45521 52216 31832

There is really just this '-997' thing in it. I think the NAs are ozeans, so is there something wrong with the data loading or do I just don't understand the data?

Best Answer

You are almost right:

NODATA is set to -32768 for oceans. Additionally, -997 is set for great lakes that are not excluded by the coastline.

Since the pixel content (growing period) makes no sense on lakes, you can safely treat -997 as NODATA too.

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Figure 2: Brazil Raster plot

Time in my laptop:

# Output:
# user  system elapsed 
# 0.203   0.033   8.688

More info about parallelization in R:

[GIS] R raster: align two rasters with different extent

You can do it by cropping the +180 -> +360 chunk out of the raster, setting its extent to -180 -> 0, and then merging it onto the 0 -> +180 chunk of the raster which you crop out with another extent:

These are the three 180-degree wide extents that we need:

> e0 = extent(raster(xmn=-180,xmx=0,ymn=-90,ymx=90))
> e1 = extent(raster(xmn=0,xmx=180,ymn=-90,ymx=90))
> e2 = extent(raster(xmn=180,xmx=360,ymn=-90,ymx=90))

Here's our raster from 0 to 360:

> extent(r2)
class       : Extent 
xmin        : 0 
xmax        : 360 
ymin        : -90 
ymax        : 90

Crop the part that is from 180 to 360, ie the farthest east (right):

> r2_e = crop(r2,e2)

Set its extent to be -180 to 0, ie on the west (left):

> extent(r2_e)=e0

Now crop out the other half of the raster, and merge it with the relocated half:

> e2_a = merge(r2_e,crop(r2,e1))

Plot and check:

> plot(e2_a)
> dim(e2_a)
[1] 360 720   1

> e2_a
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 2.683839e-06, 0.9999899  (min, max)

You should probably make sure there's no single-pixel errors around the join - it seems to have the right dimension but I've not inspected this finely.

Best Answer

Related Solutions

[GIS] Processing vector to raster faster with R

Load libraries and example data:

Simple thread example

Multithread thread example

[GIS] R raster: align two rasters with different extent

Related Question