Project

General

Profile

Download (9.96 KB) Statistics
| Branch: | Revision:
1 77ed568c Jim Regetz
library(raster)
2
3
#Get list of tile names and add full pathname to beginning
4
Tiles <- list.files("DEM/asterGdem/N59to81_W20toE19", pattern="^ASTGTM.*_dem.tif$")
5
Tiles<- paste("DEM/asterGdem/N59to81_W20toE19/", Tiles, sep="")
6
7
# datavalues: 0= water, -9999= nodata
8
9
#Create new vectors for name of tile and accompanying % of cells= ocean and non-landmass
10
#   IF % OCEAN AND % NON-LAND ARE VERY DIFFERENT USE % NON-LAND FOR QC CHECK, IF THEY
11
#   ARE THE SAME, % OCEAN IS FASTER TO CALCULATE SO USE THIS
12
13
l<- length (Tiles)
14
Names<- c(1:l)*NA
15
PctNotLand<- c(1:l) * NA
16
PctOcean<- c(1:l)*NA
17
18
#Fill vectors with % of each tiles where pixel value  = 0 (water) OR 0 & -9999 (water or nodata)
19
for (i in 1:l){
20
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
21
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
22
    PctNotLand[i]<- (count(raster(Tiles[i]),0)+count(raster(Tiles[i]),-9999))/ncell(raster(Tiles[i]))
23
}
24
25
data.frame (Names, PctOcean,PctNotLand)
26
27
#Use the above calculation to calculate % landmass for each tile and compare total over all tiles
28
#  in folder to % landmass of USGS coverage for same area.
29
30
# Requires aggregating aster tiles to match pixel size of USGS coverage:
31
32
LandAgg<- c(1:l) * NA
33
Counts_Agg<- c(1:l)*NA
34
for (i in 1:l){
35
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
36
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
37
}
38
39
a<-sum(LandAgg)
40
41
42
43
#Calculate % landmass in USGS coverage and compare:
44
45
AstRast<- raster("DEM/asterGdem/N59to81_W20toE19/w020n90/W020N90_clipped.dem")
46
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
47
1-(b/a)  #Difference= 4.49%
48
#------------------------------------------------------------------
49
#after first try, anytime PctNotLand was >= 99%, so too was PctOcean and visa versa
50
#    All subsequent runs will only calc. PctOcean to save time
51
52
Tiles <- list.files("DEM/asterGdem/N59to81_E20to59", pattern="^ASTGTM.*_dem.tif$")
53
Tiles<- paste("DEM/asterGdem/N59to81_E20to59/", Tiles, sep="")
54
55
l<- length (Tiles)
56
Names<- c(1:l)*NA
57
PctNotLand<- c(1:l) * NA
58
PctOcean<- c(1:l)*NA
59
for (i in 1:l){
60
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
61
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
62
}
63
64
data.frame (Names, PctOcean)
65
66
LandAgg<- c(1:l) * NA
67
Counts_Agg<- c(1:l)*NA
68
for (i in 1:l){
69
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
70
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
71
    print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
72
}
73
74
a<-sum(LandAgg)
75
AstRast<- raster("DEM/asterGdem/N59to81_E20to59/e020n90/E020N90_clipped.dem")
76
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
77
1-(b/a) #Difference= 2.21%
78
#-----------------------------------------------------------------------
79
Tiles <- list.files("DEM/asterGdem/N59to81_E60to99", pattern="^ASTGTM.*_dem.tif$")
80
Tiles<- paste("DEM/asterGdem/N59to81_E60to99/", Tiles, sep="")
81
82
l<- length (Tiles)
83
Names<- c(1:l)*NA
84
PctNotLand<- c(1:l) * NA
85
PctOcean<- c(1:l)*NA
86
for (i in 1:l){
87
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
88
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
89
}
90
91
data.frame (Names, PctOcean)
92
93
LandAgg<- c(1:l) * NA
94
Counts_Agg<- c(1:l)*NA
95
for (i in 1:l){
96
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
97
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
98
}
99
100
a<-sum(LandAgg)
101
AstRast<- raster("DEM/asterGdem/N59to81_E60to99/USGS_ErosDEM_N59to81E60to99/E060N90_clipped.dem")
102
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
103
1-(b/a)  # Datasets differ by 2.27%
104
105
cellStats(AstRast,max)
106
cellStats(AstRast,min)
107
count(AstRast, NA)
108
#-----------------------------------------------------------------------
109
Tiles <- list.files("DEM/asterGdem/N59to81_E100to139", pattern="^ASTGTM.*_dem.tif$")
110
Tiles<- paste("DEM/asterGdem/N59to81_E100to139/", Tiles, sep="")
111
112
l<- length (Tiles)
113
Names<- c(1:l)*NA
114
PctNotLand<- c(1:l) * NA
115
PctOcean<- c(1:l)*NA
116
for (i in 1:l){
117
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
118
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
119
}
120
121
data.frame (Names, PctOcean)
122
123
LandAgg<- c(1:l) * NA
124
Counts_Agg<- c(1:l)*NA
125
for (i in 1:l){
126
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
127
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
128
    print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
129
}
130
131
a<-sum(LandAgg)
132
AstRast<- raster("DEM/asterGdem/N59to81_E100to139/e100n90/E100N90_clipped.dem")
133
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
134
1-(b/a)  #Difference= 3.1%
135
#-----------------------------------------------------------------------
136
Tiles <- list.files("DEM/asterGdem/N59to81_E140to179", pattern="^ASTGTM.*_dem.tif$")
137
Tiles<- paste("DEM/asterGdem/N59to81_E140to179/", Tiles, sep="")
138
139
l<- length (Tiles)
140
Names<- c(1:l)*NA
141
PctNotLand<- c(1:l) * NA
142
PctOcean<- c(1:l)*NA
143
for (i in 1:l){
144
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
145
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
146
}
147
148
data.frame (Names, PctOcean)
149
150
LandAgg<- c(1:l) * NA
151
Counts_Agg<- c(1:l)*NA
152
for (i in 1:l){
153
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
154
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
155
    print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
156
}
157
158
a<-sum(LandAgg)
159
AstRast<- raster("DEM/asterGdem/N59to81_E140to179/e140n90/E140N90_clipped.dem")
160
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
161
1-(b/a)   #Difference= 3.4%
162
#-----------------------------------------------------------------------
163
Tiles <- list.files("DEM/asterGdem/N59to81_W60to21", pattern="^ASTGTM.*_dem.tif$")
164
Tiles<- paste("DEM/asterGdem/N59to81_W60to21/", Tiles, sep="")
165
166
l<- length (Tiles)
167
Names<- c(1:l)*NA
168
PctNotLand<- c(1:l) * NA
169
PctOcean<- c(1:l)*NA
170
for (i in 1:l){
171
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
172
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
173
}
174
175
data.frame (Names, PctOcean)
176
177
LandAgg<- c(1:l) * NA
178
Counts_Agg<- c(1:l)*NA
179
for (i in 1:l){
180
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
181
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
182
    print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
183
}
184
185
a<-sum(LandAgg)
186
AstRast<- raster("DEM/asterGdem/N59to81_W60to21/w060n90/W060N90_clipped.dem")
187
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
188
1-(b/a)  #Difference= 3.03%
189
#-----------------------------------------------------------------------
190
Tiles <- list.files("DEM/asterGdem/N59to81_W180to141", pattern="^ASTGTM.*_dem.tif$")
191
Tiles<- paste("DEM/asterGdem/N59to81_W180to141/", Tiles, sep="")
192
193
l<- length (Tiles)
194
Names<- c(1:l)*NA
195
PctNotLand<- c(1:l) * NA
196
PctOcean<- c(1:l)*NA
197
for (i in 1:l){
198
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
199
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
200
}
201
202
data.frame (Names, PctOcean)
203
204
LandAgg<- c(1:l) * NA
205
Counts_Agg<- c(1:l)*NA
206
for (i in 1:l){
207
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
208
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
209
     print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
210
}
211
212
a<-sum(LandAgg)
213
AstRast<- raster("DEM/asterGdem/N59to81_W180to141/w180n90/W180N90_clipped.dem")
214
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
215
1-(b/a) #Difference= 2.87%
216
#-----------------------------------------------------------------------
217
Tiles <- list.files("DEM/asterGdem/N59to81_W140to101", pattern="^ASTGTM.*_dem.tif$")
218
Tiles<- paste("DEM/asterGdem/N59to81_W140to101/", Tiles, sep="")
219
220
l<- length (Tiles)
221
Names<- c(1:l)*NA
222
PctNotLand<- c(1:l) * NA
223
PctOcean<- c(1:l)*NA
224
for (i in 1:l){
225
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
226
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
227
}
228
229
data.frame (Names, PctOcean)
230
231
LandAgg<- c(1:l) * NA
232
Counts_Agg<- c(1:l)*NA
233
for (i in 1:l){
234
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
235
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
236
    print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
237
}
238
239
a<-sum(LandAgg)
240
AstRast<- raster("DEM/asterGdem/N59to81_W140to101/w140n90/W140N90_clipped.dem")
241
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
242
1-(b/a) #Difference= 3.84%
243
#-----------------------------------------------------------------------
244
Tiles <- list.files("DEM/asterGdem/N59to81_W100to61", pattern="^ASTGTM.*_dem.tif$")
245
Tiles<- paste("DEM/asterGdem/N59to81_W100to61/", Tiles, sep="")
246
247
l<- length (Tiles)
248
Names<- c(1:l)*NA
249
PctNotLand<- c(1:l) * NA
250
PctOcean<- c(1:l)*NA
251
for (i in 1:l){
252
    Names[i]<-sapply(strsplit(Tiles[i], '[/]'), '[[', 4)
253
    PctOcean[i]<- count(raster(Tiles[i]),0)/ncell(raster(Tiles[i]))
254
}
255
256
data.frame (Names, PctOcean)
257
258
LandAgg<- c(1:l) * NA
259
Counts_Agg<- c(1:l)*NA
260
for (i in 1:l){
261
  Counts_Agg[i]<- count(aggregate(raster(Tiles[i]),30),0)+count(aggregate(raster(Tiles[i]),30),-9999)
262
  LandAgg[i]<-  ncell(aggregate(raster(Tiles[i]),30))-Counts_Agg[i]
263
    print( paste(round(100*i/l),"%",sep=""), quote=FALSE )
264
}
265
266
a<-sum(LandAgg)
267
AstRast<- raster("DEM/asterGdem/N59to81_W100to61/w100n90/W100N90_clipped.dem")
268
b<-ncell(AstRast)-cellStats(AstRast, "countNA") #Total # cells- cells with no value
269
1-(b/a) #Difference= 7.21% (lot's of coastline in this set of tiles!)
270
#-------------------------------------------------------------------
271
c<- raster("DEM/asterGdem/N59to81_E60to99/USGS_ErosDEM_N59to81E60to99/E060N90.DEM")
272
a<- raster("DEM/asterGdem/N59to81_E60to99/USGS_ErosDEM_N59to81E60to99/E060N90_clipped.dem")
273
b<- raster(nrow=82800, ncol=144000)
274
d<- raster("DEM/asterGdem/N59to81_E60to99/ASTGTM_N59E060_dem.tif")
275
s<- resample(a,b,method='ngb')
276
disaggregate(a,30)
277
e<-aggregate(d,30)
278
279
count(e, -9999)
280
count (e,0)