Project

General

Profile

Download (6.3 KB) Statistics
| Branch: | Revision:
1
# CHECK PART1
2

    
3
# R script for verifying the ASTER GDEM2 tilenames. This script will 
4
# verify that lower lefthand coordinate of tiles matches filenames
5
# (this is how the files are supposed to be named) in both dem.tif and
6
# num.tif files. 
7
# The result of Checks can be found at 
8
#    "/data/project/organisms/DEM/Yuni/documents/check/check_result.pdf" 
9
#
10
#
11
# Original author: Natalie Robinson
12
# [08-Nov-2011]
13
#    Edits by Jim and Yuni, focusing on streamlining and improving
14
#    runtime efficiency. [9-Dec-2011]
15
#    If you would like to separate the bad-name tiles in the separate 
16
#    folder, please see Nathalie's original script, located at 
17
#    "/data/project/organisms/DEM/asterGdem/R_files/AsterCheck_demAndnum.r."
18

    
19

    
20

    
21

    
22
--------------------------------------------------------------------
23
library(rgdal)
24

    
25
# path to base directory containing the tiles
26
aster.dir <- "/data/project/organisms/DEM/asterGdem2"
27

    
28
# get list of relevant ASTER GDEM tile names
29
aster.list <- list.files(aster.dir, pattern="^ASTGTM2_.*.tif$")
30

    
31

    
32

    
33

    
34
# for a given file, return TRUE if file name matches the internally
35
# specified lower left corner, otherwise print the bad-named tiles
36
# with the correct name. 
37

    
38
check.name <- function(tile_num.name, path=".", silent=TRUE) {
39
    # build expected filenameche
40
    origin <- round(GDALinfo(file.path(path, tile_num.name),
41
        silent=silent)[c("ll.x", "ll.y")])
42
    ly <- origin["ll.y"]
43
    y <- sprintf("%s%02d", if (ly>=0) "N" else "S", abs(ly))
44
    lx <- origin["ll.x"]
45
    x <- sprintf("%s%03d", if (lx>=0) "E" else "W", abs(lx))
46

    
47
    tile.suffix <- substring(tile_num.name, 16)
48
    expected.name <- paste("ASTGTM2_", y, x, tile.suffix, sep="")
49

    
50
    # compare to actual filename
51
    if (tile_num.name==expected.name) {
52
        TRUE
53
    } else {
54
        print ("Bad Name:", tile_num.name, "but supposed to be",expected.name)
55
    }
56
}
57

    
58

    
59
 
60
title_num.check <- sapply(aster.list, check.name, path=aster.dir)
61
bad.data <- data.frame(expected=title_num.check[title_num.check!="TRUE"])
62
bad.data   #this will print list of bad data 
63

    
64

    
65

    
66

    
67

    
68

    
69

    
70
-------------------------------------------------------------------
71

    
72
# CHECK PART2
73

    
74
# R script for calculating the fraction of land mass in each large USGS
75
# elevation tile, and comparing this to the fraction of land mass across
76
# all ASTER tiles covering the same spatial extent.
77
#
78
# Note that this script assumes all ASTER pixels that have a value of 0
79
# are ocean cells. This is indeed how ocean is masked out in GDEM, but
80
# we do not attempt to distinguish these from legitimate land cells that
81
# happen to have an elevation of 0 meters. We just assume the latter are
82
# rare enough to be inconsequential...
83
#
84
# Original author: Natalie Robinson
85
# [08-Nov-2011]
86
#    Edits by Jim and Yuni, focusing on streamlining and improving
87
#    runtime efficiency of the first batch of code (i.e., for
88
#    N59to81_W20toE19), which can now easily be generalized to a
89
#    function that can be reused for all other regions. [9-Dec-2011]
90
#
91
# The result of Checks can be found at 
92
#    "/data/project/organisms/DEM/Yuni/documents/check/check_result.pdf" 
93

    
94

    
95
library(raster)
96
library(rgdal) 
97

    
98
# path to base directory containing the tiles
99
aster.dir <- "/data/project/organisms/DEM/Yuni/Data/aster2"
100
usgs.dir <-  "/data/project/organisms/DEM/Yuni/Data/GTOPO30/"
101

    
102
# get list of relevant ASTER GDEM tile names and list of clipped USGS dem
103
aster.clip <- list.files(aster.dir, pattern="^aster2_.*_82N.tif$")
104
usgs.clip  <- list.files(usgs.dir) 
105

    
106

    
107

    
108
compare <- function(aster.tile.names, usgs.tile.name, aster.dir=".", usgs.dir=".") {
109

    
110
# summarize all aster tiles
111
    p.not.land <- t(sapply(aster.tile.names, function(tile.name, dir) {
112
        elev <- values(raster(file.path(dir, tile.name)))
113
        p.ocean <- mean(elev==0)
114
        p.nodata <- mean(elev==-9999)
115
        c(p.ocean=p.ocean, p.nodata=p.nodata)
116
    }, dir=aster.dir))
117

    
118
    # calculate overall proportion of all ASTER pixels containing known land
119
    # elevations (i.e., neither ocean nor nodata)
120
    aster.land.fraction <- 1 - mean(rowSums(p.not.land))
121

    
122

    
123
    # summarize usgs tiles:
124
    # calculate proportion of land in the corresponding USGS elevation
125
    # raster, in which all non-land pixels are automatically read in as NA
126
    usgs <- values(raster(file.path(usgs.dir, usgs.tile.name)))
127
    usgs.land.fraction <- mean(!is.na(usgs))
128

    
129
    # compare them:
130
    # what's the proportional difference between USGS and ASTER?
131
    diff.aster.usgs <- 1-(usgs.land.fraction/aster.land.fraction)
132

    
133
    # print the comparison details 
134
    list(
135
    comparison = c(aster.land.fraction=aster.land.fraction,
136
    usgs.land.fraction=usgs.land.fraction,
137
    difference=diff.aster.usgs),
138
    aster.details = p.not.land
139
    )
140
}
141

    
142

    
143

    
144
num.regions <- 9
145
comparisons <- lapply(1:num.regions, function(region, aster.dir, usgs.dir) {
146
    aster.tiles <- aster.clip[region]
147
    usgs.tile <- usgs.clip[region]
148
    compare(aster.tiles, usgs.tile, aster.dir, usgs.dir)
149
}, aster.dir=aster.dir, usgs.dir=usgs.dir)
150

    
151

    
152

    
153

    
154

    
155
--------------------------------------------------------------------
156

    
157
# CHECK PART3
158

    
159
# R script for calculating the number of tiles for GDEM2 and GDEM1.
160
# Quick and rough comparison of DEM tile holdings on eos for ASTER GDEM1
161
# vs ASTER GDEM2, based purely on file names
162
#
163
# Jim Regetz
164
# NCEAS
165
# Created on 09-Dec-2011
166

    
167

    
168

    
169
# list all *_dem.tif files (full paths starting from their respective
170
# base directories)
171
gdem1 <- system('find /data/project/organisms/DEM/asterGdem -name "*_dem.tif"', intern=TRUE)
172
gdem2 <- system('find /data/project/organisms/DEM/asterGdem2 -name "*_dem.tif"', intern=TRUE)
173

    
174
# extract file names, and remove any duplicates
175
# note: getting rid of spurious "._ASTGTM2xxx" files from GDEM2 list
176
f1 <- unique(basename(gdem1))
177
f2 <- unique(grep("\\._", basename(gdem2), value=TRUE, invert=TRUE))
178

    
179
# parse out latitude (always N here) and longitude values
180
a1 <- data.frame(
181
    lat = as.numeric(sub(".*N(..).*", "\\1", f1)),
182
    lon = sub(".*N..(....).*$", "\\1", f1),
183
    row.names = f1)
184
a2 <- data.frame(
185
    lat = as.numeric(sub(".*N(..).*", "\\1", f2)),
186
    lon = sub(".*N..(....).*$", "\\1", f2),
187
    row.names = f2)
188

    
189
# count files by latitude, merge results, and compute differences in
190
# counts
191
comp <- merge(data.frame(gdem1=table(a1$lat)),
192
    data.frame(gdem2=table(a2$lat)), by=1, all=TRUE)
193
comp[2:3][is.na(comp[2:3])] <- 0
194
names(comp) <- c("latitude", "gdem1count", "gdem2count")
195
comp$diff <- comp$gdem2count - comp$gdem1count
196

    
197

    
(5-5/5)