.Rhistory

setwd(oldwd)
return(NULL)
}
# Sort list of downloaded files by date
dates = extractDate(basename(files),asDate=TRUE, pos1=10, pos2=16)$inputLayerDates
observationsbydate = data.frame(dates,files)
observationsbydate = by(as.data.frame(observationsbydate), as.data.frame(observationsbydate)[,"dates"], function(x) x)
# Process all Tiles of each observation date
for (k in 1:nrow(observationsbydate)) {
filesvalid=TRUE #help variable in case a HDF File is corrupted
HDFlistbydate =  observationsbydate[[k]]
GTifflist <- vector(length=nrow(HDFlistbydate))
GTifflist2 <- vector(length=nrow(HDFlistbydate))
# Go through each Tile of the current date and: check file integrity, extract EVI layer, cut to shapefile, reproject to EPSG:4326
for (i in 1:nrow(HDFlistbydate)){
cat('Processing ... ',as.character(HDFlistbydate$files[i]),'\n',sep='')
GTifflist[i]=tempfile('NDVI',tempfolder,fileext = ".tif")
GTifflist2[i]=tempfile('NDVI',tempfolder, fileext = ".tif")
t <- tryCatch(gdalinfo(HDFlistbydate$files[i]), warning = function(w) {NULL} ,error = function(e) {NULL})
if (is.null(t)) {
unlink(as.character(HDFlistbydate$files[i]))
try(getHdf(HdfName=basename(as.character(HDFlistbydate$files[i])), wait=10))
t <- tryCatch(gdalinfo(HDFlistbydate$files[i]), warning = function(w) {NULL} ,error = function(e) {NULL})
if (is.null(t)) {
cat('Download was not successfull. Try to manually download the following file and check its validity: ', as.character(HDFlistbydate$files[i]),'\n',sep='')
filesvalid=FALSE
setwd(oldwd)
break
}
}
sds<-getSds(as.character(HDFlistbydate$files[i]))
qualityband<-raster(readGDAL(sds$SDS4gdal[12], as.is=TRUE, silent=TRUE))
mask<-calc(qualityband, fun = function(x) {bitwAnd(2,x)})
rm(qualityband);gc() # free memory
evi<-raster(readGDAL(sds$SDS4gdal[2], as.is=TRUE, silent=TRUE))
if (cloudmask) {
evi[mask==2]<- -10000
evi[is.na(evi)]<- -10000
}
rm(mask);gc()
writeRaster(evi,filename=GTifflist[i],format="GTiff",datatype="INT2S")
rm(evi);gc()
gdalwarp(srcfile=GTifflist[i],dstfile=GTifflist2[i],cutline=shapefilepath,crop_to_cutline = TRUE, t_srs="EPSG:4326") #Transform GTiff and crop to shapefile
}
# In case all Tiles of the current date are valid, mosaic them to one Gtiff File
if (filesvalid) {
filename=paste(HDFlistbydate$dates[i],'.tif',sep='')
dstfile <- file.path(dstfolder,filename)
mosaic_rasters(GTifflist2,dstfile,co=compressionmethod)  #
outputdates=c(outputdates,as.character(HDFlistbydate$dates[i]))
outputfiles=c(outputfiles, dstfile)
}
}
output=data.frame(outputfiles,outputdates)
names(output)[1] <- 'file'
names(output)[2] <- 'date'
# Delete all temporary working data and HDF Files if argument hdfstorage is NULL (no persistent storage)
do.call(unlink, list(tempfolder,recursive=TRUE))
if (is.null(hdfstorage)) {
delHdf("MOD13Q1",ask=FALSE)
delHdf("MYD13Q1",ask=FALSE)
}
setwd(oldwd)
return(output)
}
isPartOfExtent <- function(outside,inside) {
outside <- round(outside,3) # to avoid error due to precision errors
inside <- round(inside,3)
left <- xmin(inside) >= xmin(outside)
right <- xmax(inside) <= xmax(outside)
top <- ymin(inside) >= ymin(outside)
bottom <- ymax(inside) <= ymax(outside)
return(all(c(left,right,top,bottom)))
}
cropFromGeotiff <- function(date, shapefilepath, srcfolder, dstfolder, compression=FALSE) {
compressionmethod <- GenerateCompressionArgument(compression)
if (!dir.exists(dstfolder)) {
dir.create(dstfolder)
}
gtifffiles = list.files(srcfolder, pattern = "\\.tif$", full.names = TRUE)
dates=as.Date(sub(".tif","",basename(gtifffiles)))
availabledata = data.frame(gtifffiles, dates)
availabledata <- availabledata[as.numeric(availabledata$dates) %in% as.numeric(seq(from=date[1],to=date[2],by=1)),]
outputdates=c()
outputfiles=c()
if (nrow(availabledata)>0) {
# Check if bounding box of destination is within bbox of source
shp <- readOGR(shapefilepath, verbose=FALSE)
r <- raster(as.character(availabledata$gtifffiles[1]))
ispartof <- isPartOfExtent(extent(r),extent(shp))
rm(r,shp)
gc()
if (ispartof) {
for (i in 1:nrow(availabledata)) {
srcfile <- availabledata$gtifffiles[i]
dstfile <- file.path(dstfolder,paste(availabledata$dates[i],'.tif',sep=''))
cat('Processing ... ',as.character(srcfile),'\n',sep='')
gdalwarp(srcfile=srcfile,dstfile=dstfile,cutline=shapefilepath,crop_to_cutline = TRUE, t_srs="EPSG:4326",co=compressionmethod)
outputdates=c(outputdates,as.character(availabledata$dates[i]))
outputfiles=c(outputfiles, dstfile)
}
} else {
warning(paste(shapefilepath,"is not within the bounding box of",availabledata$gtifffiles[1]),sep=" ")
}
return(outputfiles)
}
output=data.frame(file=outputfiles,date=outputdates)
if (nrow(output)==0) {
return(NULL)
} else {
return(output)
}
}
check_available_data <- function(datapath, timeseries_filename) {
csv_file <- file.path(datapath,timeseries_filename)
if (file.exists(csv_file)) {
ts <- read.csv(csv_file)
ts_dates <- as.data.frame(as.Date(ts$date))
} else {
ts_dates <- NULL
}
geotiff_list <- list.files(datapath, pattern="*.tif")
if (length(geotiff_list>0)) {
geotiff_dates=as.data.frame(as.Date(sub(".tif","",basename(geotiff_list))))
} else {
geotiff_dates <- NULL
}
return(list(TS = unname(ts_dates), GEOTIFF = unname(geotiff_dates)))
}
get_latest_observation <- function(datapath, geotiff=FALSE, timeseries_filename) {
available_dates <- check_available_data(datapath, timeseries_filename)
if (geotiff) {
available_dates <- available_dates$GEOTIFF
} else {
available_dates <- available_dates$TS
}
if (is.null(available_dates)) {
return(NULL)
} else {
return(as.Date(max(unlist(available_dates))))
}
}
UpdateAndProcess <- function(database, storage_location, modis_datastorage=NULL, timeseries_filename = "timeseries.csv", max_download_chunk=15, geotiff_compression = TRUE) {
# Check Storage Location for RAW Modis Data (or set up a temporary folder) and output files and stop, if they do not exist.
if (!isString(modis_datastorage)) {
localArcPath <- file.path(tempdir(), "MODIS")
dir.create(localArcPath)
cat("Persistent Storage of MODIS RAW Data is turned off. Configurate MODIS_DATASTORAGE to a valid path to enable it.")
} else if (dir.exists(modis_datastorage)) {
localArcPath <- modis_datastorage
} else {
stop(paste("The path modis_datastorage=",modis_datastorage," does not exist",sep=""))
}
if (!dir.exists(storage_location)) {
stop(paste("The path storage_location=",storage_location," does not exist",sep=""))
}
# Initialise MODIS package.
MODISoptions(MODISserverOrder="LAADS",quiet=TRUE,localArcPath=localArcPath,outDirPath=storage_location)
# Check and Rearrange database list. Order dataframe such that subregions come last.
if (!is.logical(database$store_geotiff) || any(is.na(database$store_geotiff))) {
stop("Invalid entry in the database for store_geotiff. Only logicals (TRUE/FALSE) are allowed.")
} else if (!is.numeric(database$store_length) | !all(database$store_length>0, na.rm=TRUE)) {
stop("Invalid entry in the database for store_length. Only numerics (1,2,3,... or NA) are allowed.")
} else if (!is.logical(database$cloud_correct) || any(is.na(database$cloud_correct))) {
stop("Invalid entry in the database for cloud_correct. Only logicals (TRUE/FALSE) are allowed.")
} else if (any(duplicated(database$ID))) {
stop("The ID entries in the given database are not unique")
} else if (any(duplicated(database$name))) {
stop("The name entries in the given database are not unique")
} else if (!isDate(database$earliestdate)) {
stop("The entries for earliestdate must be of format YYYY-MM-DD or NA (NA -> as far back in time as possibe)")
} else if (!isDate(database$latestdate)) {
stop("The entries for latestdate must be of format YYYY-MM-DD or NA (NA -> today")
} else if (!all(database$is_subregion_of[isString(database$is_subregion_of)] %in% database$ID)) {
stop("One or more entries for is_subregion_of do not have a correspondend entry for a parent region")
} else if (!all(file.exists(as.character(database$shapefile)))) {
stop("One or more of the shapefiles in the database do not exist. Make sure the pathname is correct")
}
database <- database[order(database$is_subregion_of, na.last =  FALSE),]
# Resolve nested subregion dependancy to the last parentregion and check for circular dependancy
for (i in 1:nrow(database)) {
parentregion=data.frame()
subregion <- as.character(database$is_subregion_of[i])
subregionlist <- c(subregion)
while (isString(subregion)) {
parentregion <- database[database$ID==subregion,]
subregion <- as.character(parentregion$is_subregion_of)
subregionlist <- c(subregionlist,subregion)
if (any(duplicated(subregionlist))) {
stop("There is a circular dependancy of subregions!")
}
}
if (nrow(parentregion)==1) {
database$is_subregion_of[i] <- as.character(parentregion$ID)
}
}
# Find the latest observation for each database entry
# Earliest date is either the entry in the database if no geotiffs or timeseries is found in the datapath. Otherwise the latest observation
# of the timeseries is taken (store_geotiff=FALSE) or the latest observation that exists as geotiff (store_geotiff=TRUE)
# Latest date is either the entry latestdate in the database, but if it is NA, latest date is set to today (systemtime)
df_dates <- data.frame()
for (i in 1:nrow(database)) {
name=as.character(database$name[i])
datapath <- file.path(storage_location,name)
date <- get_latest_observation(datapath, geotiff=database$store_geotiff[i], timeseries_filename=timeseries_filename)
if (is.null(date)) {
startdate <- as.Date(database$earliestdate[i])
} else {
startdate <- date+1
}
datevalue <- database$latestdate[i]
if (is.na(datevalue)) {
enddate <- Sys.Date()
} else {
enddate <- as.Date(database$latestdate[i])
}
df <- data.frame(ID = database$ID[i], startdate=startdate, enddate = enddate)
df_dates <- rbind(df_dates,df)
}
rownames(df_dates) <- df_dates$ID
# daterange from earliest to latest date of all database entries
startdate <- min(df_dates$startdate)
enddate <- max(df_dates$enddate)
daterange = c(startdate,enddate)
daterange_days = enddate-startdate
# Split processing window into chunks if the daterange exceeds maxDOWNLOADchunk
if (daterange_days>max_download_chunk) {
chunks_startdate <- seq(startdate,enddate, by=max_download_chunk)
} else {
chunks_startdate <- daterange[1]
}
downloadchunks <- data.frame(start=chunks_startdate, end=c(chunks_startdate[-1]-1,enddate))
# Start Updating data: Loop over all downloadchunks resp. daterange pieces. After every chunk, delete temporary files to free harddisk space.
# Within the outer loop, loop over every entry of the database
for (j in 1:nrow(downloadchunks)) {
removefinally <- c() # Empty character vector that collects temporary files, that are not required anymore
for (i in 1:nrow(database)) {
ID <- as.character(database$ID[i])
# if startdate of database entry is within downloadchunk window, begin updating data
# crop daterange if shapefiles startdate/enddate is later/earlier than startdate/enddate of downloadchunk
if ((downloadchunks$end[j] >= df_dates[ID,"startdate"]) && (downloadchunks$start[j] <= df_dates[ID,"enddate"])) {
if (df_dates[ID,"startdate"] > downloadchunks$start[j]) {
daterange[1] <- df_dates[ID,"startdate"]
} else {
daterange[1] <- downloadchunks$start[j]
}
if (df_dates[ID,"enddate"] < downloadchunks$end[j]) {
daterange[2] <- df_dates[ID,"enddate"]
} else {
daterange[2] <- downloadchunks$end[j]
}
# concentate datapath from entry name and storage location.
# Create folder if does not yet exist
name=as.character(database$name[i])
datapath <- file.path(storage_location,name)
if (!dir.exists(datapath)) {
dir.create(datapath, recursive=TRUE)
}
# Now start downloading and processing new MODIS observations
cat('\n','############### Data for ',name,' are being updated from ',as.character(daterange[1]),' to ',as.character(daterange[2]),' ... ###############','\n',sep='')
shapefilepath <- as.character(database$shapefile[i])
# fetch the entry of the parentregion if current entry is a subregion.
subregion <- as.character(database$is_subregion_of[i])
parentregion <- database[database$ID==subregion,]
# If the current entry is a a subregion, fetch data from parent region datapath. Otherwise access MODIS FTP via MODIS package
if (nrow(parentregion)==1) {
cat('... using data from PARENTREGION with ID ',as.character(parentregion$ID),' ...','\n',sep='')
srcdatapath <- file.path(storage_location,as.character(parentregion$name))
rasterimages <- cropFromGeotiff(date = daterange, shapefilepath = shapefilepath, srcfolder = srcdatapath, dstfolder = datapath, compression = geotiff_compression)
} else {
cat('... using data from MODIS FTP server ...','\n',sep='')
rasterimages <- getMODISNDVI(date = daterange, shapefilepath=shapefilepath, dstfolder=datapath, hdfstorage=modis_datastorage, compression=geotiff_compression) #Download&Process MODIS Data
}
# Add new observations to timeseries file and tag rasterimages, that are not anymore required.
if (!is.null(rasterimages)) {
cat('Updating Time Series for ',name,' ...','\n',sep='')
# Get a list of all available rasterimages/geotiffs. Extract date vector from filenames and create dataframe with filename-date pairs.
gtifffiles <- list.files(datapath, pattern = "\\.tif$", full.names = TRUE)
dates=as.Date(sub(".tif","",basename(gtifffiles)))
datainstorage = data.frame(gtifffiles, dates)
# Read any existing timeseries file, othwerwise create empty dataframe. Extract the dates from datainstorage, for which no entry in the timeseries exists.
csvpath <- file.path(datapath,timeseries_filename)
if (file.exists(csvpath)) {
ts <- read.csv(csvpath,stringsAsFactors = FALSE, header = TRUE)
newtsdata <- datainstorage[!as.Date(datainstorage$dates) %in% as.Date(ts$date),]
} else {
ts <- data.frame()
newtsdata <- datainstorage
}
# if new data for timeseries are available, read the corresponding rasterimages and add value to the timeseries.
if (nrow(newtsdata)>0) {
values <- vector(mode="numeric", length=length(newtsdata[,1]))
dates <- vector(mode='character',length=length(newtsdata[,1]))
for (k in 1:length(values)) {
r <- raster(as.character(newtsdata[k,'gtifffiles']))*0.0001
r[r==-1]=NA
values[k]=mean(values(r), na.rm=TRUE)
dates[k]=as.character(newtsdata[k,'dates'])
}
ts <- rbind(ts,data.frame(date=dates,value=values))
ts = ts[order(ts$date, decreasing=FALSE),]
write.csv(ts,file=csvpath,row.names=FALSE)
}
# Add rasterimages that shall not be stored and are thus no longer required after the current download/daterangechunk to the vector removefinally.
# If cloud correct is activated, keep the most recent rasterimage, even if store_geotiff is set to FALSE
if (!database$store_geotiff[i]) {
if (database$cloud_correct[i]) {
datainstorage2remove <- datainstorage[!datainstorage$dates %in% max(datainstorage$dates),] #Excludes the most recent files
removefinally <- c(removefinally,as.character(datainstorage2remove$gtifffiles))
} else {
removefinally <- c(removefinally,list.files(datapath, pattern = "\\.tif$", full.names = TRUE))
}
} else {
if (!is.na(database$store_length[i]) && is.numeric(database$store_length[i]) && database$store_length[i]>0) {
alldates <- datainstorage$dates
alldates <- alldates[order(alldates,decreasing=TRUE)]
dates2keep <- alldates[1:round(database$store_length[i])]
datainstorage2remove <- datainstorage[!datainstorage$dates %in% dates2keep,]
removefinally <- c(removefinally,as.character(datainstorage2remove$gtifffiles))
}
}
} else {
cat('No new data were found for ',name,'\n',sep='')
}
}
}
if (length(removefinally)>0) {
file.remove(removefinally)
}
}
}
# Read list of shapefiles
database <- database <- read.csv(DATABASE_LOC, comment.char='#', stringsAsFactors = TRUE, colClasses = c("character","character","character","character","logical","numeric","logical","character","character"))
UpdateAndProcess(database, storage_location=DATASTORAGE_LOC, modis_datastorage = MODIS_DATASTORAGE, max_download_chunk = maxDOWNLOADchunk)
debugSource('~/Desktop/MODIS_Service/download.R')
debugSource('~/Desktop/MODIS_Service/download.R')
debugSource('~/Desktop/MODIS_Service/download.R')
MODIS::getCollection("MOD10A1")
MODIS::getCollection("MOD13Q1")
MODIS:::collections
MODIS:::collections
MODIS:::collections
MODIS:::MODIS_Products
MODIS:::MODIS_Products[1]
MODIS:::MODIS_Products[2]
MODIS:::MODIS_Products[,2]
MODIS:::MODIS_Products[2,]
MODIS:::MODIS_Products
MODIS::getProduct("MOD13Q1")
MODIS:::addProduct(product="MCD64A1", sensor = "MODIS", platform = "Combined", pf1="MOTA",
pf2="MOD", topic="Burned Area", type = "Tile", res ="500m", temp_res="Monthly",
internalseparator = "\.", server ="LPDAAC",path_ext = "/Users/gogunkoya/Downloads/MODIS/external", overwrite = TRUE)
MODIS:::MODIS_Products
install.packages("MODISSnow")
library(MODISSnow)
download_data
download_data()
download_data("2017-01-01")
download_data(as.Date("2017-01-01"))
download_data(as.Date("2017-01-01"),h = 18, v = 5)
download_data(as.Date("2017-01-01"),h = 18, v = 5, user="test", password="test")
MODIS:::MODIS_FTPinfo
getProduct("MOD16")
getProduct("MOD10")
getProduct()
download_data(as.Date("2017-01-01"),h = 18, v = 5, user="julesair2", password="538-FuJnS")
download_data(as.Date("2017-01-01"),h = 18, v = 5, user="julesair2", password="538-FuJnS")
MODIS:::addServer(name = "NSIDC", sensor = "MODIS",
basepath = "https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/2018.02.20/",
varpath = "YYYY.MM.DD/")
MODIS:::MODIS_FTPinfo
MODIS:::addServer(name = "NSIDC", sensor = "MODIS",
basepath = "https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/2018.02.20/",
varpath = "PF1/PRODUCT.CCC/DATE/")
/home/jules/R/x86_64-pc-linux-gnu-library/3.4/MODIS/external
MODIS:::addServer(name = "NSIDC", sensor = "MODIS",
basepath = "https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/2018.02.20/",
varpath = "PF1/PRODUCT.CCC/DATE/", path_ext = "/home/jules/R/x86_64-pc-linux-gnu-library/3.4/MODIS/external")
MODIS:::MODIS_FTPinfo
MODIS:::addServer(name = "NSIDC", sensor = "MODIS",
basepath = "https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/2018.02.20/",
varpath = "PF1/PRODUCT.CCC/DATE/", path_ext = "/home/jules/R/x86_64-pc-linux-gnu-library/3.4/MODIS/external", overwrite=TRUE)
MODIS:::MODIS_FTPinfo
load("/home/jules/R/x86_64-pc-linux-gnu-library/3.4/MODIS/external/MODIS_FTPinfo.RData")
library(MODIS)
MODIS:::MODIS_FTPinfo
library(RCurl)
filesUrl
filesUrl()
library(rvest)
install.packages("rvest")
library(rvest)
read_html(https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/)
read_html("https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/")
getURI("https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/")
library(RCurl)
library(XML)
XML::getHTMLLinks("https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/")
getURI("https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/")
write("test.html",getURI("https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/"))
write(getURI("https://n5eil01u.ecs.nsidc.org/MOST/MOD10A1.006/"),"test.html")
getURI("https://n5eil01u.ecs.nsidc.org/opendap/MOST/MOD10A1.006/2000.02.24/contents.html")
install.packages("httr")
install.packages("httr")
library(hrrt)
library(httr)
GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS"))
)
GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS"))
write("test.html",GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS")))
out <- GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS"))
write(out,"test.html")
out
1
write(out,"test2.html")
write(out$content,"test2.html")
GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS"),  write_disk(destfile="test.html")
)
GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS"), write_disk("test.html"))
GET("https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/",
authenticate("julesair2", "538-FuJnS"), write_disk("test.html",overwrite=TRUE))
readLines(out)
str(content(out))
out$content
str(out$content)
content(out, "text")
htmlout <- content(out, "text")
links <- xpathSApply(htmlout, "//a/@href")
library(XML)
test<-xmlTreeParse(htmlout)
test<-xmlTreeParse("test.html")
test<-htmlTreeParse("test.html")
test
library(xml2)
library(rvest)
head(html_attr(html_nodes(htmlout, "a"), "href"))
head(html_attr(html_nodes(out, "a"), "href"))
read("test.html")
readline("test.html")
readline("test.html")
readLines("test.html")
htmlout <- content(out, "text")
page=read_html("http://www.yelp.com/search?find_loc=New+York,+NY,+USA")
htmlParse(out)
links <- xpathSApply(htmlParse(out), "//a/@href")
links
links[1]
links[2]
links[3]
links[4]
links[5]
links[6]
links[7]
links[8]
links[7]
links[7]
as.Date(links[7])
as.Date(links[7],format=yyyy.mm.dd/)
as.Date(links[7],format="yyyy.mm.dd/")
as.Date(links[7],format="%Y.%m.%d/")
as.Date(links[1],format="%Y.%m.%d/")
date <- as.Date(links[7],format="%Y.%m.%d/")
nextlink <- paste(baselink,links[7])
baselink <- "https://n5eil01u.ecs.nsidc.org/MOSA/MYD10A1.005/"
nextlink <- paste(baselink,links[7])
nextlink
nextlink <- paste(baselink,links[7],sep="")
nextlink
out <- GET(nextlink,
authenticate("julesair2", "538-FuJnS"))
links2 <- xpathSApply(htmlParse(out), "//a/@href")
links2[1]
links2[2]
links2[3]
links2[4]
links2[5]
links2[6]
links2[7]
links2[8]
links2[9]
links2[1000]
links2[1001]
source('~/Desktop/extract_nsdic.R')
source('~/Desktop/extract_nsdic.R')
source('~/Desktop/extract_nsdic.R')
debugSource('~/Desktop/MODIS_Service/download.R')
debugSource('~/Desktop/MODIS_Service/download.R')
hist(r)
hist(r)
hist(r)
test <- raster(imagefile)
image(test[test==-32766])
image(test[test==-32767])
image(test==-32766)
image(old_r)
image(old_r==-32767)
image(old_r==-32766)
image(old_r==-32766)
image(old_r==-32767)
image(old_r==-32765)
image(r)
image(r==-32767)
image(r==-32766)
image(r==-32767)
image(old_r==-32767)
image(r==-32767)