License New BSD license
Lines 148
Keywords
extents (3) GDAL (2) GIS (7) OGR (2) overlap (1)
Permissions
Viewable by Everyone
Editable by All Siafoo Users
Hide
Siafoo – the intersection of pastebin, help desk, version control, and social networking Join Siafoo Now or Learn More

Find GIS Files that Overlap Extents / Find Extents of GIS File Atom Feed 0

In Brief Returns a list of GIS files that overlap the given extents, or finds the extents of any GDAL- or OGR-compatible GIS data.... more
# 's
  1#!/usr/bin/env python
2
3# Copyright 2007 Regents of the University of California
4# Written by David Isaacson at the University of California, Davis
5# BSD License
6
7import gdal,ogr,getopt,sys,os
8from gdalconst import *
9
10def main():
11 gdal.PushErrorHandler()
12
13 #define some defaults
14 path = ""
15 suffixes = ""
16 delimiter = " "
17 minx = -180
18 miny = -90
19 maxx = 180
20 maxy = 90
21
22 try:
23 opts,args = getopt.getopt(sys.argv[1:],"d:s:")
24 except getopt.GetoptError:
25 #print help information and exit:
26 usage()
27 sys.exit(2)
28
29 for opt,arg in opts:
30 if opt == "-s":
31 suffixes = arg.split()
32 if opt == "-d":
33 delimiter = arg
34
35 if len(args) == 5:
36 path = os.path.abspath(args[0])
37 minx = long(args[1])
38 miny = long(args[2])
39 maxx = long(args[3])
40 maxy = long(args[4])
41 elif len(args) == 1:
42 path = os.path.abspath(args[0])
43 else:
44 usage()
45 sys.exit(2)
46
47 if os.path.isdir(path):
48 overlaps = findOverlaps(path,suffixes,minx,miny,maxx,maxy)
49 print delimiter.join(overlaps)
50 elif os.path.isfile(path):
51 coords = findCoordinates(path)
52 if coords:
53 print "%.8Lf %.8Lf %.8Lf %.8Lf" %tuple(coords)
54 else:
55 print "The path",path,"does not seem to exist, or is not a regular file or directory."
56 sys.exit();
57
58def findOverlaps(path,suffixes,sminx,sminy,smaxx,smaxy):
59 minx = long(sminx);
60 miny = long(sminy);
61 maxx = long(smaxx);
62 maxy = long(smaxy);
63 overlaps = []
64
65 rasters = getFilesBySuffix(path,suffixes)
66 for raster in rasters:
67 coords = findCoordinates(raster)
68 if not coords:
69 continue
70 if extentsOverlap([minx,miny,maxx,maxy],coords):
71 overlaps.append(raster);
72 return overlaps
73
74def findOverlapsFromVector(path,suffixes,vector):
75 return findOverlaps(path,suffixes,vector[0],vector[1],vector[2],vector[3])
76
77def findCoordinates(path):
78 coords = findGDALCoordinates(path)
79 if not coords:
80 coords = findOGRCoordinates(path)
81 if not coords:
82 return []
83 return coords
84
85def findGDALCoordinates(path):
86 if not os.path.isfile(path):
87 return []
88 data = gdal.Open(path,GA_ReadOnly)
89 if data is None:
90 return []
91 geoTransform = data.GetGeoTransform()
92 minx = geoTransform[0]
93 maxy = geoTransform[3]
94 maxx = minx + geoTransform[1]*data.RasterXSize
95 miny = maxy + geoTransform[5]*data.RasterYSize
96 return [minx,miny,maxx,maxy]
97
98def findOGRCoordinates(path):
99 if not os.path.isfile(path):
100 return[]
101 dataSource = ogr.Open(path)
102 if not dataSource:
103 return []
104 extents = []
105 for i in range(0,dataSource.GetLayerCount()):
106 layer = dataSource.GetLayer(i)
107 extents.append(layer.GetExtent())
108
109 dataSource.Release()
110
111 if not extents:
112 return []
113 extent = extents[0]
114 minx = extent[0]
115 maxx = extent[1]
116 miny = extent[2]
117 maxy = extent[3]
118
119 if len(extents) > 1:
120 for extent in extents[1:]:
121 if extent[0]<minx: minx = extent[0]
122 if extent[1]>maxx: maxx = extent[1]
123 if extent[2]<miny: miny = extent[2]
124 if extent[3]>maxy: maxy = extent[3]
125
126 return [minx,miny,maxx,maxy]
127
128
129def extentsOverlap(a,b):
130 if a[0]<=b[2] and b[0]<=a[2] and a[1]<=b[3] and b[1]<=a[3]:
131 return 1
132 else:
133 return 0
134
135def getFilesBySuffix(path,suffixes):
136 rasters = []
137 if os.path.isdir(path):
138 os.path.walk(path,checkSuffixes,[rasters,suffixes])
139 elif os.path.isfile(path):
140 for suffix in suffixes:
141 if path.endswith("." + suffix):
142 rasters.append(path)
143 break
144 return rasters
145
146def checkSuffixes(array,path,files):
147 rasters,suffixes = array;
148 for file in files:
149 if suffixes:
150 for suffix in suffixes:
151 suffix = "." + suffix;
152 if file.endswith(suffix):
153 rasters.append(os.path.join(path,file))
154 break
155 else:
156 rasters.append(os.path.join(path,file))
157
158def usage():
159 print "Usage:"
160 print " myraster.py [-d \"delimiter\"] [-s \"suffix1 suffix2\"] directory_name [xmin ymin xmax ymax]"
161 print " myraster.py file_name"
162
163if __name__ == "__main__":
164 main()

Returns a list of GIS files that overlap the given extents, or finds the extents of any GDAL- or OGR-compatible GIS data.

Incorporates snippets Get Files With Suffixes, Find Extents of GIS Raster Data, and Find Extents of GIS Vector Data