Subset HDF5 matrix using values from another submatrix

Question

I have been working with HDF5 files and I have been able to do some subsetting in rhdf5. there are three files:. Longitude, latitude and the ColumnAmountNO2Trop. I have extracted this for all the days in the year as a list in "files".

files <- list.files(pattern = ".he5", full.names = TRUE)
attribute <- "/HDFEOS/SWATHS/ColumnAmountNO2/Data Fields/ColumnAmountNO2Trop"
attribute2<-"/HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Longitude"
attribute3<-"/HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Latitude"

and the sub files below:

out.list <- lapply(files, h5read, attribute)
Lon <- lapply(files, h5read, attribute2)
Lat<-lapply(files, h5read, attribute3)

However, I need to subset out.list(which contains the 'ColumnAmountNO2Trop' for all the days of the year) based on the latitude and longitude values to narrow down my area of reference geographically. I was able to subset them using row and column numbers:

lapply(out.list, function(x) x[2:8,2:8]) 

However, the goegraphic location of 2,2 on day one may not be the same on day 2. I tried to define the Longitude and latitude values to subset with below but it returned an error message.

Lond<-c(2,9)
Latd<-c(2,9)
lonKeep <- which(Lon > Lond[1] & Lon < Lond[2])
latKeep <- which(lat> latRan[1] & lat< latRan[2]) 

How do I subset 'out.list' for Lon 2-9 and Lat 2-9, please?

Answer 1

You'll find other approaches, possibly with a solution closer to you needs. An option with hdf5 and raster would be to extract the relevant data from hdf5 files, build a raster, crop it to the ROI and get the values for that area.

I'd do something like this:

library(raster)
library(maptools)
source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")
library(rhdf5)
library(latticeExtra)

my_wd <- './Stackoverflow/22474417'
files <- list.files(path = my_wd, pattern = ".he5", full.names = F)
files
#[1] "M1.he5" "M2.he5"

attribute <- "/HDFEOS/SWATHS/ColumnAmountNO2/Data Fields/ColumnAmountNO2Trop"
attribute2<- "/HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Longitude"
attribute3<- "/HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Latitude"

Read a single file

m1 <- h5read(file.path(my_wd, files[1]), name = attribute)
dim(m1) # file dimension
# [1] 60 54
prod(dim(m1))
# [1] 3060

We'll use it to build a rasterLayer, after extracting the geographical extent from atribute2 and atribute3

Lon <- h5read(file.path(my_wd, files[1]) , attribute2)
Lat <- h5read(file.path(my_wd, files[1]) , attribute3)

xmin <- min(Lon[1:prod(dim(m1))]) # Min. Longitude
# [1] -7.141283
xmax <- max(Lon[1:prod(dim(m1))]) # Max. Longitude
ymin <- min(Lat[1:prod(dim(m1))]) # Min. Longitude
ymax <- max(Lat[1:prod(dim(m1))]) # Max. Longitude

We can build a raster with the info above

m1m <- matrix(m1, nrow = 60)    
m1r <- raster(m1m, xmn = xmin, xmx = xmax,
              ymn =  ymin, ymx = ymax)

Take some spatial data to overlay

data(wrld_simpl)
spdata <- wrld_simpl[which(wrld_simpl@data$NAME %in% c('Nigeria', 'Cameroon', 'Benin',
                                                       'Togo', 'Ghana',"Cote d'Ivoire",
                                                       'Gabon', 'Equatorial Guinea')), ] 

From Africa Shoreline 30m

delta <- readOGR(dsn = './africa_shoreline_30m',
                 layer = 'nigeria_delta')

Build a ROI extent

frm <- extent(c(2, 9, 2, 9))
pfrm <- as(frm, 'SpatialPolygons')

plot it

spplot(m1r,scales = list(draw = TRUE),  ylim=c(-1, 10)) +
  latticeExtra::layer(sp.polygons(stp, fill = NA, col = 'blue'))+
  latticeExtra::layer(sp.polygons(pfrm, fill = NA, col = 'red'))

Crop and get values from ROI

m1rf <- crop(m1r, frm)

spplot(m1rf, scales = list(draw = TRUE), xlim = c(1, 10), ylim=c(1, 10)) +
  latticeExtra::layer(sp.lines(delta, fill = NA, col = 'blue'))+
  latticeExtra::layer(sp.polygons(pfrm, fill = NA, col = 'red'))

summary(m1rf)
                layer
Min.    -6.528723e+15
1st Qu.  9.437798e+14
Median   1.440395e+15
3rd Qu.  1.896734e+15
Max.     4.232078e+15
NA's     0.000000e+00

m1vals <- getValues(m1rf)

Once you agree with this, it's easy to loop over your file folder and get your data.