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
Stay up to dateembedded code automagically updates, each snippet and article has a feed 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