Sorry but I'm pretty new in RStudio, I couldn't even do reprex But I'm sending the whole code below and erased some unnecessary parts. I'm trying to create a date list by using POSIXct function, it should be an interval between 2015-04-01 and 2015-04-30 but couldn't become successful. In this original code, I think user tried to get rid of some specific days such as 14, 17 etc. In my own code, I don't need this but the code itself is not working even in its original form.
# Make_SMAP_L3_36km_SM_raster_stack
library(gdalUtils)
library(rgdal)
library(raster)
library(matlab)
library(stringr)
wd <- "C:/SM_code/"
filenames<- "D:/SMAP_L3_SM_P_E_20150425_R16020_001.h5" # Need to be modified to your directory
library(rhdf5)
h5ls(filenames)
SMAP.sm_am <- h5read(filenames, "Soil_Moisture_Retrieval_Data_AM/soil_moisture")
SMAP.SM <- h5read(filenames, "Soil_Moisture_Retrieval_Data_PM/soil_moisture_pm")
SMAP.lat <- h5read(filenames, "Soil_Moisture_Retrieval_Data_PM/latitude_pm")
SMAP.lon <- h5read(filenames, "Soil_Moisture_Retrieval_Data_PM/longitude_pm") #soil_moisture: 1- 483, soil_moisture_pm: 484 -
II<-which(SMAP.lat == -9999.00000)
SMAP.lat[II] <- NaN
SMAP.lat1 <- SMAP.lat
SMAP.lat1[,] <- NaN
for (i in 1:1624/4){
I <- which(is.na(SMAP.lat[,i])==FALSE)
SMAP.lat1[,i] <- SMAP.lat[I[1],i]}
SMAP.lon1 <- SMAP.lon
SMAP.lon1[,] <- NaN
II<-which(SMAP.lon1 == -9999.00000)
SMAP.lon1[II] <- NaN
for (i in 1:3856/4){
I <- which(is.na(SMAP.lon[i,])==FALSE)
SMAP.lon1[i,] <- SMAP.lon[i,I[1]]}
library(sp)
library(rgdal)
nlat<-47
lat_len<-34
wlon<-201
lon_len<-62
slat<-80
elon<-262
# NCRFC extent: -104.8587, -82.42284, 37.49958, 50.31208 (xmin, xmax, ymin, ymax)
# indices: 201,262,80,47 (length 62,34)
SMAP.lon2<-SMAP.lon1[wlon:elon,nlat:slat]
SMAP.lat2<-SMAP.lat1[wlon:elon,nlat:slat]
LON <- unlist(as.data.frame(SMAP.lon2))
LAT <- unlist(as.data.frame(SMAP.lat2))
coords1 <- as.data.frame(cbind(LON,LAT)) #make data frame object with coordinate vectors
points <- SpatialPoints(coords1) #make spatial object with lat/long coordinates
#set projection of coordinates to lat/long
proj4string(points)<-CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#crs(points)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#tranform lat/lon to EASE-Grid 2.0 (original EASE, global cylindrical equal area, WGS84 projection)
points<-spTransform(points,CRS("+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30 +ellps=WGS84 +datum=WGS84 +units=m"))
# Read SMAP SM
yb<-paste("D:/SMAP_Data", sep="")
SMAP.list<-list.files(yb)
SMAP.list<-unlist(str_extract_all(SMAP.list,".+(.h5)"))
# Make date list
SMAP_datetime <- as.POSIXct(paste0(substr(SMAP.list,14,17),"-",substr(SMAP.list,18,19),"-",substr(SMAP.list,20,21)))
# Select specific periods e.g. 2015-03-01 ~ 2015-05-31
startdate <- "2015-04-12"
enddate <- "2015-04-20" # 03-27
start <- which(SMAP_datetime == startdate)
end <- which(SMAP_datetime == enddate)
N <- end - start + 1; N
# check actual period
as.POSIXct(enddate) - as.POSIXct(startdate)
# Expand RRB boundary
rfc_us_ease@bbox[1] <- extent(rfc_us_ease)@xmin - 50000
rfc_us_ease@bbox[3] <- extent(rfc_us_ease)@xmax + 50000
rfc_us_ease@bbox[2] <- extent(rfc_us_ease)@ymin -50000
rfc_us_ease@bbox[4] <- extent(rfc_us_ease)@ymax +50000
#SM<-repmat(NA,lon_len,lat_len, N) # SMAP.list
for (ampm in 1:2) {
if (ampm == 1) {overpass1 <- "AM"; overpass <- "am"}
if (ampm == 2) {overpass1 <- "PM"; overpass <- "pm"}
for (i in start:end) { #1:length(SMAP.list)
j <- i-start+1
if (ampm == 1) {SMAP.SM <- h5read(paste0(yb, SMAP.list[i]), paste0("Soil_Moisture_Retrieval_Data_AM/soil_moisture"))
sm_error <- h5read(paste0(yb, SMAP.list[i]), paste0("Soil_Moisture_Retrieval_Data_AM/soil_moisture_error"))} ######### AM / PM
if (ampm == 2) {SMAP.SM <- h5read(paste0(yb, SMAP.list[i]), paste0("Soil_Moisture_Retrieval_Data_PM/soil_moisture_pm"))
sm_error <- h5read(paste0(yb, SMAP.list[i]), paste0("Soil_Moisture_Retrieval_Data_PM/soil_moisture_error_pm"))} ######### AM / PM
######### AM / PM
II<-which(SMAP.SM == -9999.00000)
SMAP.SM[II] <- NaN
#SM[,,j]<-SMAP.SM[wlon:elon,nlat:slat]
SMSM <-SMAP.SM[wlon:elon,nlat:slat]
LL<-which(sm_error == -9999.00000)
sm_error[LL] <- NaN
err <-sm_error[wlon:elon,nlat:slat]
data<-as.data.frame(err)
data<-unlist(data) #unlist data frame to get vector of soil moisture
data <- SpatialPointsDataFrame(points, data.frame(data)) #add soil moisture data to spatial object
###instead, used "SpatialPixelsDataFrame" with tolerance limit for irregularity
###see bottom of: http://gis.stackexchange.com/questions/79062/how-to-make-raster-from-irregular-point-data-without-interpolation
data <- SpatialPixelsDataFrame(data, tolerance = 8.17194e-05, data@data) # 8.17194e-05
data <- raster(data[,'data'])
if(i==start) {
SMAP_L3_36km_raster<-data
} else {SMAP_L3_36km_raster<-addLayer(SMAP_L3_36km_raster,data)} #create stack of raster files
print(j)
SMAP_datetime <- seq(as.Date(startdate),as.Date(enddate), by=1)
if (ampm == 1 && i == end) {
SMAP_L3_36km_AM_USraster <- SMAP_L3_36km_raster
save(SMAP_L3_36km_AM_USraster, SMAP_datetime, file = paste0("D:/Ongoing/SMAP_L3_36km_US_",overpass,"_",startdate,"_",enddate,".RData"))
}
if (ampm == 2 && i == end) {SMAP_L3_36km_PM_USraster <- SMAP_L3_9km_raster
save(SMAP_L3_36km_PM_USraster, SMAP_datetime, file = paste0("D:/Ongoing/SMAP_L3_36km_US_",overpass,"_",startdate,"_",enddate,".RData"))
}
}
}
library(rasterVis)
levelplot(SMAP_L3_36km_PM_USraster, layers = 1:3, par.settings = RdBuTheme, main = paste0("SMAP"))#+ layer(sp.polygons(rfc_us_ease))