-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathDownload_NAVI.R
140 lines (101 loc) · 3.71 KB
/
Download_NAVI.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
library(cptcity)
library(raster)
library(stars)
library(rgee)
library(sf)
library(rgdal)
library(tiff)
library(googledrive)
library(stringr)
ee_Initialize(drive = TRUE)
#ee_check()
#Define a region of interest
roi <- ee$Geometry$Polygon(
proj = "EPSG:4326",
geodesic = FALSE,
coords = list(
c(19.5696, -33.660),
c(18.146, -34.418),
c(19.5696, -34.418),
c(18.146, -33.660)
)
)
modis_ndvi <- ee$ImageCollection("MODIS/006/MOD13A2")
#MODIS makes it simple to filter out poor quality pixels thanks to a quality control bits band (DetailedQA).
#The following function helps us to distinct between good data (bit == …00) and marginal data (bit != …00).
getQABits <- function(image, qa) {
# Convert binary (character) to decimal (little endian)
qa <- sum(2^(which(rev(unlist(strsplit(as.character(qa), "")) == 1))-1))
# Return a mask band image, giving the qa value.
image$bitwiseAnd(qa)$lt(1)
}
#Using getQABits we construct a single-argument function (mod13A2_clean) that is used to map over all the images of the collection (modis_ndvi).
mod13A2_clean <- function(img) {
# Extract the NDVI band
ndvi_values <- img$select("NDVI")
# Extract the quality band
ndvi_qa <- img$select("SummaryQA")
# Select pixels to mask
quality_mask <- getQABits(ndvi_qa, "11")
# Mask pixels with value zero.
ndvi_values$updateMask(quality_mask)
}
#Filter the collection (modis_ndvi) by a date range.
ndvi_composite <- modis_ndvi$
filter(ee$Filter$date('2001-01-01', '2001-02-01'))$
#filter(ee$Filter$lte("CLOUDY_PIXEL_PERCENTAGE", 20))$
map(mod13A2_clean)
#OPTIONAL: Use Map to display the results in an interactive way.
#It only works for display one image, still needs to develop
#Do not run it yet
scale <- 0.0001
Map$setCenter(lon =18.83606, lat = -34.06176, zoom = 7)
Map$addLayer(
eeObject = ndvi_composite,
visParams = list(
min = 0.2 / scale,
max = 0.7 / scale,
palette = cpt("grass_ndvi", 10)
)
) + Map$addLayer(roi)
# checking the exsisting datasets in google drive
# path --- a single folder on Google Drive whose contents you want to list
# pattern --- items whose names match this regular expression are returned.
raster_find <- drive_ls(pattern = "20",path='rgee_backup')
# collecting the date of exisiting datasets
existing_date <- as.character(as.list(raster_find$name))
existing_date <- substr(existing_date,1,10)
# collecting the date of datasets you want to export
desire_date <- as.list(ee_get_date_ic(ndvi_composite)[,1])
desire_date <- substr(desire_date,19,29)
# a loop to download data
`%!in%` <- Negate(`%in%`)
for (i in 1:length(desire_date)){
if (desire_date[i] %!in% existing_date){
missing_date<-desire_date[i] # data that you want to export is existing on the folder
missing_date <- str_replace_all(missing_date,'_','-')
print(missing_date)
ndvi_composite <- modis_ndvi$
filter(ee$Filter$date( missing_date,as.character(as.Date(missing_date)+1)))$ # Specifies the date
map(mod13A2_clean)
mod_ndvi <- ee_imagecollection_to_local(
ic = ndvi_composite,
scale = 100,
region = roi,
via = 'drive'
)
}
}
# access to your google drive
# path --- a single folder on Google Drive whose contents you want to list
# pattern --- items whose names match this regular expression are returned.
raster_find <- drive_ls(pattern = "2001_01_17",path='rgee_backup')
raster_path <- drive_get(raster_find$name)
# download data to your directory
# setwd()
# getwd()
raster_download<-drive_download(raster_path$name,type = "tif",overwrite = TRUE)
# produce raster objects
imported_raster <- raster(raster_download$local_path, as.is=TRUE)
imported_raster <- aggregate(imported_raster,fact=2)
plot(imported_raster)