| 1 | from osgeo import * |
|---|
| 2 | import osgeo.ogr |
|---|
| 3 | import osgeo.gdal |
|---|
| 4 | import libxml2 |
|---|
| 5 | import os |
|---|
| 6 | import sys |
|---|
| 7 | import zoo |
|---|
| 8 | |
|---|
| 9 | def readFileFromBuffer(data,ext): |
|---|
| 10 | try: |
|---|
| 11 | geometry=[] |
|---|
| 12 | print >> sys.stderr,'/vsimem//temp1'+ext |
|---|
| 13 | #print >> sys.stderr,data |
|---|
| 14 | osgeo.gdal.FileFromMemBuffer('/vsimem//temp1'+ext,data) |
|---|
| 15 | ds = osgeo.ogr.Open('/vsimem//temp1'+ext) |
|---|
| 16 | lyr = ds.GetLayer(0) |
|---|
| 17 | feat = lyr.GetNextFeature() |
|---|
| 18 | while feat is not None: |
|---|
| 19 | geometry+=[feat.Clone()] |
|---|
| 20 | feat.Destroy() |
|---|
| 21 | feat = lyr.GetNextFeature() |
|---|
| 22 | ds.Destroy() |
|---|
| 23 | osgeo.gdal.Unlink('/vsimem//temp1'+ext) |
|---|
| 24 | return geometry |
|---|
| 25 | except Exception,e: |
|---|
| 26 | print >> sys.stderr,e |
|---|
| 27 | return [] |
|---|
| 28 | |
|---|
| 29 | def buildFeatureFromGeomtry(conf,geom,driverName,ext): |
|---|
| 30 | drv = osgeo.ogr.GetDriverByName( driverName ) |
|---|
| 31 | ds = drv.CreateDataSource( "/vsimem//store"+conf["lenv"]["sid"]+"0."+ext ) |
|---|
| 32 | lyr = ds.CreateLayer( "Result", None, osgeo.ogr.wkbUnknown ) |
|---|
| 33 | field_defn = osgeo.ogr.FieldDefn( "Name", osgeo.ogr.OFTString ) |
|---|
| 34 | field_defn.SetWidth( len("Result10000") ) |
|---|
| 35 | lyr.CreateField ( field_defn ) |
|---|
| 36 | feat = osgeo.ogr.Feature(lyr.GetLayerDefn()) |
|---|
| 37 | feat.SetField( "Name", "Input0" ) |
|---|
| 38 | feat.SetGeometry(geom) |
|---|
| 39 | lyr.CreateFeature(feat) |
|---|
| 40 | ds.Destroy() |
|---|
| 41 | return [feat] |
|---|
| 42 | |
|---|
| 43 | def createGeometryFromWFS(conf,my_wfs_response): |
|---|
| 44 | try: |
|---|
| 45 | geom=osgeo.ogr.CreateGeometryFromGML(my_wfs_response.replace('<?xml version="1.0" encoding="utf-8"?>\n','')) |
|---|
| 46 | except: |
|---|
| 47 | geom=None |
|---|
| 48 | try: |
|---|
| 49 | if geom is None: |
|---|
| 50 | if not(conf["lenv"].has_key("cnt")): |
|---|
| 51 | conf["lenv"]["cnt"]=0 |
|---|
| 52 | else: |
|---|
| 53 | conf["lenv"]["cnt"]+=1 |
|---|
| 54 | return readFileFromBuffer(my_wfs_response,str(conf["lenv"]["cnt"])) |
|---|
| 55 | else: |
|---|
| 56 | return buildFeatureFromGeomtry(conf,geom,"GML","xml") |
|---|
| 57 | except: |
|---|
| 58 | print >> sys.stderr,"Unable to load file input data !!!\n\n\n" |
|---|
| 59 | |
|---|
| 60 | def createLayerFromJson(conf,obj): |
|---|
| 61 | geom=osgeo.ogr.CreateGeometryFromJson(obj) |
|---|
| 62 | if geom is None: |
|---|
| 63 | return readFileFromBuffer(obj,".json") |
|---|
| 64 | else: |
|---|
| 65 | return buildFeatureFromGeomtry(conf,geom,"GeoJSON","json") |
|---|
| 66 | |
|---|
| 67 | def extractInputs(conf,obj): |
|---|
| 68 | if obj["mimeType"]=="application/json": |
|---|
| 69 | return createLayerFromJson(conf,obj["value"]) |
|---|
| 70 | else: |
|---|
| 71 | return createGeometryFromWFS(conf,obj["value"]) |
|---|
| 72 | |
|---|
| 73 | def outputResult(conf,obj,geom): |
|---|
| 74 | driverName = "GML" |
|---|
| 75 | extension = [ ".xml" , ".xsd" ] |
|---|
| 76 | if obj["mimeType"]=="application/json": |
|---|
| 77 | driverName = "GeoJSON" |
|---|
| 78 | extension = [ ".js" ] |
|---|
| 79 | if obj.keys().count("schema")>0 and \ |
|---|
| 80 | obj["schema"]=="http://schemas.opengis.net/kml/2.2.0/ogckml22.xsd": |
|---|
| 81 | driverName = "KML" |
|---|
| 82 | extension = [ ".kml" ] |
|---|
| 83 | drv = osgeo.ogr.GetDriverByName( driverName ) |
|---|
| 84 | # Create virtual file |
|---|
| 85 | ds = drv.CreateDataSource( "/vsimem/store"+conf["lenv"]["sid"]+extension[0] ) |
|---|
| 86 | lyr = ds.CreateLayer( "Result", None, osgeo.ogr.wkbUnknown ) |
|---|
| 87 | i=0 |
|---|
| 88 | while i < len(geom): |
|---|
| 89 | if i==0 and driverName!="GeoJSON": |
|---|
| 90 | poDstFDefn=geom[i].GetDefnRef() |
|---|
| 91 | if poDstFDefn is not None: |
|---|
| 92 | nDstFieldCount = poDstFDefn.GetFieldCount() |
|---|
| 93 | for iField in range(nDstFieldCount): |
|---|
| 94 | poSrcFieldDefn = poDstFDefn.GetFieldDefn(iField) |
|---|
| 95 | oFieldDefn = osgeo.ogr.FieldDefn(poSrcFieldDefn.GetNameRef(),poSrcFieldDefn.GetType()) |
|---|
| 96 | oFieldDefn.SetWidth( poSrcFieldDefn.GetWidth() ) |
|---|
| 97 | oFieldDefn.SetPrecision( poSrcFieldDefn.GetPrecision() ) |
|---|
| 98 | lyr.CreateField( oFieldDefn ) |
|---|
| 99 | lyr.CreateFeature(geom[i]) |
|---|
| 100 | geom[i].Destroy() |
|---|
| 101 | i+=1 |
|---|
| 102 | ds.Destroy() |
|---|
| 103 | vsiFile=osgeo.gdal.VSIFOpenL("/vsimem/store"+conf["lenv"]["sid"]+extension[0],"r") |
|---|
| 104 | i=0 |
|---|
| 105 | while osgeo.gdal.VSIFSeekL(vsiFile,0,os.SEEK_END)>0: |
|---|
| 106 | i+=1 |
|---|
| 107 | fileSize=osgeo.gdal.VSIFTellL(vsiFile) |
|---|
| 108 | osgeo.gdal.VSIFSeekL(vsiFile,0,os.SEEK_SET) |
|---|
| 109 | obj["value"]=osgeo.gdal.VSIFReadL(fileSize,1,vsiFile) |
|---|
| 110 | osgeo.gdal.Unlink("/vsimem/store"+conf["lenv"]["sid"]+extension[0]) |
|---|
| 111 | |
|---|
| 112 | def BufferPy(conf,inputs,outputs): |
|---|
| 113 | print >> sys.stderr, "Starting service ..." |
|---|
| 114 | try: |
|---|
| 115 | bdist=float(inputs["BufferDistance"]["value"]) |
|---|
| 116 | except: |
|---|
| 117 | bdist=1 |
|---|
| 118 | print >> sys.stderr, bdist |
|---|
| 119 | geometry=extractInputs(conf,inputs["InputPolygon"]) |
|---|
| 120 | i=0 |
|---|
| 121 | rgeometries=[] |
|---|
| 122 | while i < len(geometry): |
|---|
| 123 | tmp=geometry[i].Clone() |
|---|
| 124 | resg=geometry[i].GetGeometryRef().Buffer(bdist) |
|---|
| 125 | tmp.SetGeometryDirectly(resg) |
|---|
| 126 | rgeometries+=[tmp] |
|---|
| 127 | geometry[i].Destroy() |
|---|
| 128 | resg.thisown=False |
|---|
| 129 | tmp.thisown=False |
|---|
| 130 | i+=1 |
|---|
| 131 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 132 | i=0 |
|---|
| 133 | return zoo.SERVICE_SUCCEEDED |
|---|
| 134 | |
|---|
| 135 | def BoundaryPy(conf,inputs,outputs): |
|---|
| 136 | geometry=extractInputs(conf,inputs["InputPolygon"]) |
|---|
| 137 | i=0 |
|---|
| 138 | rgeometries=[] |
|---|
| 139 | while i < len(geometry): |
|---|
| 140 | tmp=geometry[i].Clone() |
|---|
| 141 | resg=geometry[i].GetGeometryRef() |
|---|
| 142 | resg=resg.GetBoundary() |
|---|
| 143 | tmp.SetGeometryDirectly(resg) |
|---|
| 144 | rgeometries+=[tmp] |
|---|
| 145 | geometry[i].Destroy() |
|---|
| 146 | i+=1 |
|---|
| 147 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 148 | return zoo.SERVICE_SUCCEEDED |
|---|
| 149 | |
|---|
| 150 | def CentroidPy(conf,inputs,outputs): |
|---|
| 151 | geometry=extractInputs(conf,inputs["InputPolygon"]) |
|---|
| 152 | i=0 |
|---|
| 153 | rgeometries=[] |
|---|
| 154 | while i < len(geometry): |
|---|
| 155 | tmp=geometry[i].Clone() |
|---|
| 156 | resg=geometry[i].GetGeometryRef() |
|---|
| 157 | if resg.GetGeometryType()!=3: |
|---|
| 158 | resg=resg.ConvexHull() |
|---|
| 159 | resg=resg.Centroid() |
|---|
| 160 | tmp.SetGeometryDirectly(resg) |
|---|
| 161 | rgeometries+=[tmp] |
|---|
| 162 | geometry[i].Destroy() |
|---|
| 163 | i+=1 |
|---|
| 164 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 165 | return zoo.SERVICE_SUCCEEDED |
|---|
| 166 | |
|---|
| 167 | def ConvexHullPy(conf,inputs,outputs): |
|---|
| 168 | geometry=extractInputs(conf,inputs["InputPolygon"]) |
|---|
| 169 | i=0 |
|---|
| 170 | rgeometries=[] |
|---|
| 171 | while i < len(geometry): |
|---|
| 172 | tmp=geometry[i].Clone() |
|---|
| 173 | resg=geometry[i].GetGeometryRef().ConvexHull() |
|---|
| 174 | tmp.SetGeometryDirectly(resg) |
|---|
| 175 | rgeometries+=[tmp] |
|---|
| 176 | geometry[i].Destroy() |
|---|
| 177 | i+=1 |
|---|
| 178 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 179 | return zoo.SERVICE_SUCCEEDED |
|---|
| 180 | |
|---|
| 181 | |
|---|
| 182 | |
|---|
| 183 | def EnvelopePy(conf,inputs,outputs): |
|---|
| 184 | print >> sys.stderr, inputs |
|---|
| 185 | try: |
|---|
| 186 | bdist=float(inputs["BufferDistance"]["value"]) |
|---|
| 187 | except: |
|---|
| 188 | bdist=10 |
|---|
| 189 | geometry=extractInputs(conf,inputs["InputPolygon"]) |
|---|
| 190 | tmp=geometry[0].GetGeometryRef().GetEnvelope() |
|---|
| 191 | outputs["Result"]["value"]=str(tmp[0])+','+str(tmp[2])+','+str(tmp[1])+','+str(tmp[3])+','+'urn:ogc:def:crs:OGC:1.3:CRS84' |
|---|
| 192 | print >> sys.stderr,outputs["Result"] |
|---|
| 193 | return 3 |
|---|
| 194 | |
|---|
| 195 | def UnionPy(conf,inputs,outputs): |
|---|
| 196 | geometry1=extractInputs(conf,inputs["InputEntity1"]) |
|---|
| 197 | geometry2=extractInputs(conf,inputs["InputEntity2"]) |
|---|
| 198 | rgeometries=[] |
|---|
| 199 | i=0 |
|---|
| 200 | while i < len(geometry1): |
|---|
| 201 | j=0 |
|---|
| 202 | while j < len(geometry2): |
|---|
| 203 | tmp=geometry1[i].Clone() |
|---|
| 204 | resg=geometry2[j].GetGeometryRef() |
|---|
| 205 | resg=resg.Union(geometry1[i].GetGeometryRef()) |
|---|
| 206 | tmp.SetGeometryDirectly(resg) |
|---|
| 207 | if not(resg.IsEmpty()): |
|---|
| 208 | rgeometries+=[tmp] |
|---|
| 209 | j+=1 |
|---|
| 210 | geometry1[i].Destroy() |
|---|
| 211 | i+=1 |
|---|
| 212 | i=0 |
|---|
| 213 | while i < len(geometry2): |
|---|
| 214 | geometry2[i].Destroy() |
|---|
| 215 | i+=1 |
|---|
| 216 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 217 | return 3 |
|---|
| 218 | |
|---|
| 219 | def IntersectionPy(conf,inputs,outputs): |
|---|
| 220 | |
|---|
| 221 | geometry1=extractInputs(conf,inputs["InputEntity1"]) |
|---|
| 222 | geometry2=extractInputs(conf,inputs["InputEntity2"]) |
|---|
| 223 | |
|---|
| 224 | print >> sys.stderr,str(len(geometry1))+" "+str(len(geometry2)) |
|---|
| 225 | |
|---|
| 226 | rgeometries=[] |
|---|
| 227 | fids=[] |
|---|
| 228 | i=0 |
|---|
| 229 | while i < len(geometry1): |
|---|
| 230 | j=0 |
|---|
| 231 | while j < len(geometry2): |
|---|
| 232 | tmp=geometry2[j].Clone() |
|---|
| 233 | resg=geometry2[j].GetGeometryRef() |
|---|
| 234 | #resg=resg.Intersection(geometry1[i].GetGeometryRef()) |
|---|
| 235 | resg=geometry1[i].GetGeometryRef().Intersection(resg) |
|---|
| 236 | tmp.SetGeometryDirectly(resg) |
|---|
| 237 | if resg is not None and not(resg.IsEmpty()) and fids.count(tmp.GetFID())==0: |
|---|
| 238 | rgeometries+=[tmp] |
|---|
| 239 | fids+=[tmp.GetFID()] |
|---|
| 240 | else: |
|---|
| 241 | tmp.Destroy() |
|---|
| 242 | j+=1 |
|---|
| 243 | geometry1[i].Destroy() |
|---|
| 244 | i+=1 |
|---|
| 245 | i=0 |
|---|
| 246 | while i < len(geometry2): |
|---|
| 247 | geometry2[i].Destroy() |
|---|
| 248 | i+=1 |
|---|
| 249 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 250 | print >> sys.stderr,"/outputResult" |
|---|
| 251 | return 3 |
|---|
| 252 | |
|---|
| 253 | def DifferencePy(conf,inputs,outputs): |
|---|
| 254 | geometry1=extractInputs(conf,inputs["InputEntity1"]) |
|---|
| 255 | geometry2=extractInputs(conf,inputs["InputEntity2"]) |
|---|
| 256 | rgeometries=[] |
|---|
| 257 | i=0 |
|---|
| 258 | while i < len(geometry1): |
|---|
| 259 | j=0 |
|---|
| 260 | while j < len(geometry2): |
|---|
| 261 | tmp=geometry2[j].Clone() |
|---|
| 262 | resg=geometry1[i].GetGeometryRef() |
|---|
| 263 | resg=resg.Difference(geometry2[i].GetGeometryRef()) |
|---|
| 264 | tmp.SetGeometryDirectly(resg) |
|---|
| 265 | if not(resg.IsEmpty()): |
|---|
| 266 | rgeometries+=[tmp] |
|---|
| 267 | j+=1 |
|---|
| 268 | geometry1[i].Destroy() |
|---|
| 269 | i+=1 |
|---|
| 270 | i=0 |
|---|
| 271 | while i < len(geometry2): |
|---|
| 272 | geometry2[i].Destroy() |
|---|
| 273 | i+=1 |
|---|
| 274 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 275 | return 3 |
|---|
| 276 | |
|---|
| 277 | def SymDifferencePy(conf,inputs,outputs): |
|---|
| 278 | geometry1=extractInputs(conf,inputs["InputEntity1"]) |
|---|
| 279 | geometry2=extractInputs(conf,inputs["InputEntity2"]) |
|---|
| 280 | rgeometries=[] |
|---|
| 281 | i=0 |
|---|
| 282 | while i < len(geometry1): |
|---|
| 283 | j=0 |
|---|
| 284 | while j < len(geometry2): |
|---|
| 285 | tmp=geometry2[j].Clone() |
|---|
| 286 | resg=geometry1[i].GetGeometryRef() |
|---|
| 287 | resg=resg.SymmetricDifference(geometry2[i].GetGeometryRef()) |
|---|
| 288 | tmp.SetGeometryDirectly(resg) |
|---|
| 289 | rgeometries+=[tmp] |
|---|
| 290 | j+=1 |
|---|
| 291 | geometry1[i].Destroy() |
|---|
| 292 | i+=1 |
|---|
| 293 | i=0 |
|---|
| 294 | while i < len(geometry2): |
|---|
| 295 | geometry2[i].Destroy() |
|---|
| 296 | i+=1 |
|---|
| 297 | outputResult(conf,outputs["Result"],rgeometries) |
|---|
| 298 | return 3 |
|---|
| 299 | |
|---|