I am trying to convert my h5 daily files into a raster format. I converted into raster format. When I extracted my area of interest. I observed that issue is related to extent. Kindly anyone guides me on how to solve this issue. The R code and image is present in link with the attachment. Thanks
https://drive.google.com/drive/folders/18-hj2LEYWBN-uIDDTdqZ-x-WUxpCJu7H?usp=sharing
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("rhdf5")
# browseVignettes("rhdf5")
library(rhdf5)
library(sp)
library(raster)
# h5ls("reconstruction_Sierra_2015.h5")
h5ls("reconstruction_indus_CY2001.h5")
# h5readAttributes(file = "reconstruction_Sierra_2015.h5", name = "Grid")
h5readAttributes(file = "reconstruction_indus_CY2001.h5", name = "Grid")
# h5f = H5Fopen("reconstruction_Sierra_2015.h5")
h5f = H5Fopen("reconstruction_indus_CY2001.h5")
# h5f
# h5f&'Grid'
#system.time( swe <- h5f$Grid$swe )
system.time( melt <- h5f$Grid$melt )
locations <- data.frame(
lon=c(74.86764,73.48753, 74.87066 , 73.37798 , 78.82102 ,75.85160 ,75.78263 , 78.46446 ),
lat = c(35.16700, 36.25674, 36.49362, 35.21188, 34.20916, 34.48459, 35.76965, 33.23380)
)
coordinates(locations) <- ~lon+lat
proj4string(locations) <- CRS("+proj=longlat")
i = 1
a = seq(1,365,1)
for (i in 1:length(a)) {
swe180 <- swe[,,i]
b <- swe180 == 65535
# table(b)
swe180[b] <- -1
b <- swe180 > 200
# table(b)
swe180[b] <- 200
b <- swe180 < 0
# table(b)
swe180[b] <- 20
# image(swe180)
# image(swe180)
# str(swe180)
# h5readAttributes(file = "reconstruction_Sierra_2016.h5", name = "Grid")$ReferencingMatrix
RM <- h5readAttributes(file = "reconstruction_indus_CY2001.h5", name = "Grid")$ReferencingMatrix
GT <- GridTopology(c(RM[3,1], RM[3,2]+RM[1,2]*dim(swe)[1]), c(RM[2,1], -RM[1,2]), c(dim(swe)[2],dim(swe)[1]))
#GT <- GridTopology(c(RM[3,1], RM[3,2]+RM[1,2]*dim(melt)[1]), c(RM[2,1], -RM[1,2]), c(dim(melt)[2],dim(melt)[1]))
# GT <- GridTopology(c(-1.088854e+07, 4718608.3619-463.3127*1978), c(463.3127, 463.3127), c(2171,1978))
# GT
SG = SpatialGrid(GT)
# str(SG)
# proj4string(SG) <- CRS("+proj=sinu")
# str(SG)
proj4string(SG) <- CRS("+proj=utm +zone=43 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
locations_aea <- spTransform(locations, CRS(proj4string(SG)))
SGDF = SpatialGridDataFrame(SG, data.frame(melt = as.numeric(t(swe180))))
gridded(SGDF)<- TRUE
r = raster(SGDF)
plot(SGDF, axes=T)
writeRaster(r,paste0('/shared/MODIS/shastaH5SWEinR/2001_swe/',i,'.tif'))
}
## Open Raster Files and Extract Area of Interest
files <- mixedsort(list.files(path = '/shared/MODIS/shastaH5SWEinR/2001_melt/2002/', pattern='.tif$',
all.files=TRUE, full.names=TRUE,recursive = TRUE))
shp= readOGR("Hunza.shp")
e = extent(shp)
shp = spTransform(shp,CRS("+proj=utm +zone=43 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
i = 1
y <- NULL
c = as.data.frame(matrix(NA,nrow = 1, ncol = 1))
for (i in 1:length(files)) {
r1 = raster(files[i])
crs(r1) = "+proj=utm +zone=43 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(r1)
r1_mask = raster::mask(r1,shp)
plot(r1_mask,axes = TRUE,ext = extent(shp))
# Extracting Values as Data Frame
r1_extract = raster::extract(r1,shp, df=TRUE,na.rm = TRUE)
# Stroing as Raster
writeRaster(r1_mask,paste0('/shared/MODIS/shastaH5SWEinR/2001_swe/Hunza/',i,'.tif'))
c = cbind(r1_extract,y)
c1=t(c)
write.csv(c1,file = 'Hunza_SWE_2001.csv')
}