Project

General

Profile

Download (1.69 KB) Statistics
| Branch: | Revision:
1 7526fb1c Jim Regetz
#This script was used to check that the coordinates at the top left corner of each mosaiced Aster tile match those that are expected based on tile name (quick QC check that mosaicing ran correctly)
2
3
#Natalie Robinson 
4
#Created on Jan. 28, 2012 
5
6
import os
7
import osgeo
8
from osgeo import gdal
9
from osgeo.gdal import *
10
import subprocess
11
import re
12
13
path= "/data/project/organisms/DEM/asterGdem2/90m_NoPixelOffset/Mosaiced"
14
#file= path+ "N60to65E000_005.tif"
15
#ds=gdal.Open(file, GA_ReadOnly)
16
17
FileList=[]
18
Filex=[]
19
Filey=[]
20
x_min=[]
21
y_max=[]
22
23
#Simpler command if no subdirectories (use for path + "/N59to60")
24
#for dirname, dirnames, filenames in os.walk(path):
25
#    for file in filenames:
26
#      ext= os.path.splitext(file)[1]
27
#      if ext == ".tif":
28
29
#If subdirectories, list subdirectories to ignore
30
ignore='N59to60'
31
for root, dirs, files in os.walk(path):
32
   if(ignore in dirs):
33
      dirs.remove(ignore)
34
   for file in files:
35
       FileList.append(path + "/" + file)
36
       Filex.append(file[7:11])
37
       Filey.append(file[5:7])
38
       for i in range (0,len(Filex)):
39
          if re.search("W",Filex[i][0]):
40
             Filex[i]= Filex[i].replace("W","-")
41
          else:
42
             Filex[i]=Filex[i].replace("E","")
43
            
44
45
for j in range (0,len(Filex)):
46
   Filex[j]=int(Filex[j])
47
   Filey[j]=int(Filey[j])
48
49
for f in FileList:
50
       ds=gdal.Open(f)
51
       geotrans= ds.GetGeoTransform()
52
       x_min.append(geotrans[0])
53
       y_max.append(geotrans[3])
54
55
for i in range (0,len(Filex)):
56
    if Filex[i]==x_min[i]:
57
        print Filex[i], "OK"
58
    else:
59
        print Filex[i], x_min
60
61
for i in range (0,len(Filey)):
62
    if Filey[i]==y_max[i]:
63
        print Filey[i], "OK"
64
    else:
65
        print Filey[i], y_max[i]
66