1
|
#This script was used to:
|
2
|
# 1) resample Aster GDEM2 tiles into 90m resolution and shift grids to get rid of #1/2 pixel offset (so that there will be no overlapping pixels between adjacent #tiles, which were found to have different values for many overlapping pixels)
|
3
|
# 2) resample SRTM tiles to get rid of 1/2 pixel offset
|
4
|
|
5
|
#Natalie Robinson
|
6
|
#Created on Jan. 18, 2012
|
7
|
|
8
|
import os
|
9
|
import sys
|
10
|
import re
|
11
|
import string
|
12
|
import osgeo
|
13
|
from osgeo import gdal
|
14
|
from osgeo.gdalconst import *
|
15
|
import subprocess
|
16
|
from subprocess import *
|
17
|
|
18
|
#AsterGdem2------------------------------------------------------------------
|
19
|
sourceDir= "/data/project/organisms/DEM/asterGdem2"
|
20
|
outDir= "/data/project/organisms/DEM/asterGdem2/90m_NoPixelOffset"
|
21
|
|
22
|
#separate out .dem files and create lists of file names (input, output,
|
23
|
#and without fullpath for easy indexing)
|
24
|
for x in os.listdir(sourceDir):
|
25
|
if re.search("\_dem.tif$",x):
|
26
|
files.append(x)
|
27
|
|
28
|
for y in range (0,len(files)):
|
29
|
inName.append (sourceDir + "/" + files[y])
|
30
|
outName.append (outDir + "/" + files[y][8:15] + "_ReSample.tif")
|
31
|
|
32
|
|
33
|
#Create lists of x_min, x_max, y_min, y_max values for automating gdalwarp:
|
34
|
files = []
|
35
|
inName=[]
|
36
|
outName=[]
|
37
|
lat=[]
|
38
|
lon=[]
|
39
|
|
40
|
#x_min and y_min (lat and lon):
|
41
|
for i in range (0,len(files)):
|
42
|
lat.append (files[i][9:11])
|
43
|
lon.append (files[i][11:15])
|
44
|
for j in range (0,len(lon)):
|
45
|
if re.search("W",lon[i][0]):
|
46
|
lon[i]= lon[i].replace("W","-")
|
47
|
else:
|
48
|
lon[i]=lon[i].replace("E","")
|
49
|
|
50
|
#x_max and y_max (lat+ 1 and lon+1, but where elements are
|
51
|
#string format (for correct read into gdalwarp)
|
52
|
lat_PlusOne=[]
|
53
|
lon_PlusOne=[]
|
54
|
|
55
|
for y in range (0,len(lon)):
|
56
|
lon[y]=int(lon[y])
|
57
|
lat[y]=int(lat[y])
|
58
|
lon_PlusOne.append(lon[y]+1)
|
59
|
lat_PlusOne.append(lat[y]+1)
|
60
|
#Convert back to string
|
61
|
lon_PlusOne=["%s" % el for el in lon_PlusOne]
|
62
|
lat_PlusOne=["%s" % el for el in lat_PlusOne]
|
63
|
lon=["%s" % el for el in lon]
|
64
|
lat=["%s" % el for el in lat]
|
65
|
|
66
|
#Create new files using names from code above
|
67
|
for i in range (0,len(files)):
|
68
|
subprocess.call (["gdalwarp","-te",lon[i],lat[i],lon_PlusOne[i],lat_PlusOne[i],"-ts","1200","1200","-srcnodata","-9999","-dstnodata","-9999","-r","bilinear",inName[i], outName[i]])
|
69
|
|
70
|
#SRTM 90m---------------------------------------------------------------------
|
71
|
sourceDir= "/data/project/organisms/DEM/cgiarSrtm/SRTM_90m_ASCII_4_1"
|
72
|
outDir= "/data/project/organisms/DEM/cgiarSrtm/SRTM_90m_ASCII_4_1/Tiles_Resampled"
|
73
|
|
74
|
#separate out .asc files and create lists of file names and convert
|
75
|
#naming scheme into lat/lon info for file creation and naming (input,
|
76
|
#output, and without fullpath for easy indexing)
|
77
|
files = []
|
78
|
inName=[]
|
79
|
outName=[]
|
80
|
s=[]
|
81
|
t=[]
|
82
|
|
83
|
for x in os.listdir(sourceDir):
|
84
|
if re.search("\.asc$",x):
|
85
|
files.append(x)
|
86
|
|
87
|
#Convert SRTM naming scheme to correct lat/lon info
|
88
|
for i in range(len(files)):
|
89
|
s.append(files[i][5:7])
|
90
|
t.append(files[i][8:10])
|
91
|
s[i]=int(s[i])
|
92
|
t[i]=int(t[i])
|
93
|
|
94
|
#Create lists of x_min, x_max, y_min, y_max values for automating gdalwarp:
|
95
|
lat=[]
|
96
|
lon=[]
|
97
|
lon_PlusFive=[]
|
98
|
lat_PlusFive=[]
|
99
|
|
100
|
#append lon (y_min) w/ correct longitude based on filename, repeat for lat (x_min)
|
101
|
#Seed lon with correct starting longitude, and lat with correct starting latitude
|
102
|
files[0] # is 03_08, or lat 20, lon -170
|
103
|
lon.append(-170)
|
104
|
lat.append(20)
|
105
|
|
106
|
for k in range (0,len(s)):
|
107
|
lon.append(lon[k]+ 5 * (s[k+1]-s[k]))
|
108
|
lat.append(lat[k]- 5 * (t[k+1]-t[k]))
|
109
|
|
110
|
#x_max and y_max (lat+ 5 and lon+5, but where elements are
|
111
|
#string format (for correct read into gdalwarp)
|
112
|
for y in range (0,len(lon)):
|
113
|
lon[y]=int(lon[y])
|
114
|
lat[y]=int(lat[y])
|
115
|
lon_PlusFive.append(lon[y]+5)
|
116
|
lat_PlusFive.append(lat[y]+5)
|
117
|
#Convert back to string
|
118
|
lon_PlusFive=["%s" % el for el in lon_PlusFive]
|
119
|
lat_PlusFive=["%s" % el for el in lat_PlusFive]
|
120
|
lon=["%s" % el for el in lon]
|
121
|
lat=["%s" % el for el in lat]
|
122
|
|
123
|
#Create input and output filenames
|
124
|
for y in range (0,len(files)):
|
125
|
inName.append (sourceDir + "/" + files[y])
|
126
|
outName.append (outDir + "/" + lat[y] + "Lat" + lon[y] + "Lon" + "_ReSample.tif")
|
127
|
|
128
|
#Create new files
|
129
|
for i in range (0,len(files)):
|
130
|
subprocess.call (["gdalwarp","-te",lon[i],lat[i],lon_PlusFive[i],lat_PlusFive[i],"-ts","6000","6000","-ot", "Int16", "-srcnodata","-9999","-dstnodata","-9999","-r","bilinear",inName[i], outName[i]])
|