source: trunk/zoo-project/zoo-services/ogr/base-vect-ops-py/cgi-env/ogr_sp.py @ 940

Last change on this file since 940 was 917, checked in by djay, 6 years ago

Merge prototype-v0 branch in trunk
  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-python
File size: 13.5 KB
Line 
1# -*- coding: utf-8 -*-
2#
3# Author : Gérald FENOY
4#
5# Copyright 2009-2013 GeoLabs SARL. All rights reserved.
6#
7# Permission is hereby granted, free of charge, to any person obtaining a
8# copy of this software and associated documentation files (the
9# "Software"), to deal in the Software without restriction, including with
10# out limitation the rights to use, copy, modify, merge, publish,
11# distribute, sublicense, and/or sell copies of the Software, and to
12# permit persons to whom the Software is furnished to do so, subject to
13# the following conditions:
14#
15# The above copyright notice and this permission notice shall be included
16# in all copies or substantial portions of the Software.
17#
18# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
19# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
21# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
22# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25#
26
27from osgeo import *
28import osgeo.ogr
29import osgeo.gdal
30import libxml2
31import os
32import sys
33import zoo
34
35def readFileFromBuffer(data,ext):
36    try:
37        geometry=[]
38        print >> sys.stderr,'/vsimem//temp1'+ext
39        #print >> sys.stderr,data
40        osgeo.gdal.FileFromMemBuffer('/vsimem//temp1'+ext,data)
41        ds = osgeo.ogr.Open('/vsimem//temp1'+ext)
42        lyr = ds.GetLayer(0)
43        feat = lyr.GetNextFeature()
44        while feat is not None:
45            geometry+=[feat.Clone()]
46            feat.Destroy()
47            feat = lyr.GetNextFeature()
48        ds.Destroy()
49        osgeo.gdal.Unlink('/vsimem//temp1'+ext)
50        return geometry
51    except Exception,e:
52        print >> sys.stderr,e
53        return []
54   
55def buildFeatureFromGeomtry(conf,geom,driverName,ext):
56    drv = osgeo.ogr.GetDriverByName( driverName )
57    ds = drv.CreateDataSource( "/vsimem//store"+conf["lenv"]["sid"]+"0."+ext )
58    lyr = ds.CreateLayer( "Result", None, osgeo.ogr.wkbUnknown )
59    field_defn = osgeo.ogr.FieldDefn( "Name", osgeo.ogr.OFTString )
60    field_defn.SetWidth( len("Result10000") )
61    lyr.CreateField ( field_defn )
62    feat = osgeo.ogr.Feature(lyr.GetLayerDefn())
63    feat.SetField( "Name", "Input0" )
64    feat.SetGeometry(geom)
65    lyr.CreateFeature(feat)
66    ds.Destroy()
67    return [feat]
68
69def createGeometryFromWFS(conf,my_wfs_response):
70    try:
71        geom=osgeo.ogr.CreateGeometryFromGML(my_wfs_response.replace('<?xml version="1.0" encoding="utf-8"?>\n','').replace('<?xml version=\'1.0\' encoding="utf-8"?>\n',''))
72    except Exception,e:
73        print >> sys.stderr,"**"
74        print >> sys.stderr,e
75        geom=None
76    try:
77        print >> sys.stderr,geom is None
78        if geom is None:
79            if not(conf["lenv"].has_key("cnt")):
80                conf["lenv"]["cnt"]=0
81            else:
82                conf["lenv"]["cnt"]+=1
83            return readFileFromBuffer(my_wfs_response,str(conf["lenv"]["cnt"]))
84        else:
85            return buildFeatureFromGeomtry(conf,geom,"GML","xml")
86    except:
87        print >> sys.stderr,"Unable to load file input data !!!\n\n\n"
88
89def createLayerFromJson(conf,obj):
90    geom=osgeo.ogr.CreateGeometryFromJson(obj)
91    if geom is None:
92        return readFileFromBuffer(obj,".json")
93    else:
94        return buildFeatureFromGeomtry(conf,geom,"GeoJSON","json")
95
96def extractInputs(conf,obj):
97    if obj.keys().count("cache_file"):
98        print >> sys.stderr,obj
99        geometry=[]
100        ds = osgeo.ogr.Open(obj["cache_file"])
101        lyr = ds.GetLayer(0)
102        feat = lyr.GetNextFeature()
103        while feat is not None:
104            geometry+=[feat.Clone()]
105            feat.Destroy()
106            feat = lyr.GetNextFeature()
107        ds.Destroy()
108        return geometry   
109    if obj["mimeType"]=="application/json":
110        return createLayerFromJson(conf,obj["value"])
111    else:
112        return createGeometryFromWFS(conf,obj["value"])
113   
114def outputResult(conf,obj,geom):
115    if obj["mimeType"].count("text/xml")>0:
116        driverName = "GML"
117        extension = [ ".xml" , ".xsd" ]
118        format_list = { "2.": 'GML2', "3.1.1": 'GML3', "3.1": 'GML3Deegree', "3.2": 'GML3.2' }
119        opts=['FORMAT=%s,GML3_LONGSRS=YES',format_list["3.2"]]
120        for i in format_list:
121            if obj["mimeType"].count(i)>0:
122                opts=['FORMAT=%s,GML3_LONGSRS=YES',format_list[i]]
123    if obj["mimeType"]=="application/json":
124        driverName = "GeoJSON"
125        extension = [ ".js" ]
126        opts=None
127    if obj.keys().count("schema")>0 and \
128            obj["schema"]=="http://schemas.opengis.net/kml/2.2.0/ogckml22.xsd":
129        driverName = "KML"
130        extension = [ ".kml" ]
131        opts=None
132    drv = osgeo.ogr.GetDriverByName( driverName )
133    print >> sys.stderr,drv
134    # Create virtual file
135    ds = drv.CreateDataSource( "/vsimem/store"+conf["lenv"]["sid"]+extension[0],options = opts)
136    print >> sys.stderr,ds
137    lyr = ds.CreateLayer( "Result", None, osgeo.ogr.wkbUnknown )
138    print >> sys.stderr,lyr
139    i=0
140    print >> sys.stderr,driverName
141    print >> sys.stderr,extension
142    while i < len(geom):
143        if i==0 and driverName!="GeoJSON":
144            poDstFDefn=geom[i].GetDefnRef()
145            if poDstFDefn is not None:
146                nDstFieldCount = poDstFDefn.GetFieldCount()
147                for iField in range(nDstFieldCount):
148                    poSrcFieldDefn = poDstFDefn.GetFieldDefn(iField)
149                    oFieldDefn = osgeo.ogr.FieldDefn(poSrcFieldDefn.GetNameRef(),poSrcFieldDefn.GetType())
150                    oFieldDefn.SetWidth( poSrcFieldDefn.GetWidth() )
151                    oFieldDefn.SetPrecision( poSrcFieldDefn.GetPrecision() )
152                    lyr.CreateField( oFieldDefn )
153        try:
154            lyr.CreateFeature(geom[i])
155        except:
156            pass
157        geom[i].Destroy()
158        i+=1
159    ds.Destroy()
160    vsiFile=osgeo.gdal.VSIFOpenL("/vsimem/store"+conf["lenv"]["sid"]+extension[0],"r")
161    i=0
162    while osgeo.gdal.VSIFSeekL(vsiFile,0,os.SEEK_END)>0:
163        i+=1
164    fileSize=osgeo.gdal.VSIFTellL(vsiFile)
165    osgeo.gdal.VSIFSeekL(vsiFile,0,os.SEEK_SET)
166    obj["value"]=osgeo.gdal.VSIFReadL(fileSize,1,vsiFile)
167    osgeo.gdal.Unlink("/vsimem/store"+conf["lenv"]["sid"]+extension[0])
168
169def BufferPy(conf,inputs,outputs):
170    print >> sys.stderr, "Starting service ..."
171    try:
172        bdist=float(inputs["BufferDistance"]["value"])
173    except:
174        bdist=1
175    print >> sys.stderr, bdist
176    geometry=extractInputs(conf,inputs["InputPolygon"])
177    i=0
178    rgeometries=[]
179    while i < len(geometry):
180        tmp=geometry[i].Clone()
181        resg=geometry[i].GetGeometryRef().Buffer(bdist)
182        tmp.SetGeometryDirectly(resg)
183        rgeometries+=[tmp]
184        geometry[i].Destroy()
185        resg.thisown=False 
186        tmp.thisown=False
187        i+=1
188    outputResult(conf,outputs["Result"],rgeometries)
189    i=0
190    return zoo.SERVICE_SUCCEEDED
191
192def Clean(conf,inputs,outputs):
193    from shapely.wkb import loads
194    print >> sys.stderr, "Starting service ..."
195    features=extractInputs(conf,inputs["InputData"])
196    i=0
197    rgeometries=[]
198    while i < len(features):
199        tmp=features[i].Clone()
200        resg=features[i].GetGeometryRef()
201        if resg is not None:
202            geom = loads(resg.ExportToWkb())
203            if geom.is_valid:
204                tmp.SetGeometryDirectly(resg)
205                rgeometries+=[tmp]
206                print >> sys.stderr,"valid !"
207            else:
208                print >> sys.stderr,"invalid !"
209                print >> sys.stderr,geom.wkt[0:50]
210        features[i].Destroy()
211        resg.thisown=False 
212        tmp.thisown=False
213        i+=1
214    outputResult(conf,outputs["Result"],rgeometries)
215    i=0
216    print >> sys.stderr,"Return"
217    return zoo.SERVICE_SUCCEEDED
218
219def TransformService(conf,inputs,outputs):
220    from osgeo import osr
221    geometry=extractInputs(conf,inputs["InputData"])
222    sourceRef = osr.SpatialReference()
223    tmp=inputs["SourceCRS"]["value"].split(":")
224    sourceRef.ImportFromEPSG(int(tmp[len(tmp)-1]))
225    targetRef = osr.SpatialReference()   
226    tmp=inputs["TargetCRS"]["value"].split(":")
227    targetRef.ImportFromEPSG(int(tmp[len(tmp)-1]))
228    transform = osr.CoordinateTransformation(sourceRef, targetRef)
229    i=0
230    rgeometries=[]
231    while i < len(geometry):
232        tmp=geometry[i].Clone()
233        resg=geometry[i].GetGeometryRef()
234        resg.Transform(transform)
235        tmp.SetGeometryDirectly(resg.Clone())
236        rgeometries+=[tmp]
237        geometry[i].Destroy()
238        i+=1
239    outputResult(conf,outputs["TransformedData"],rgeometries)
240    return zoo.SERVICE_SUCCEEDED
241
242def BoundaryPy(conf,inputs,outputs):
243    geometry=extractInputs(conf,inputs["InputPolygon"])
244    i=0
245    rgeometries=[]
246    while i < len(geometry):
247        tmp=geometry[i].Clone()
248        resg=geometry[i].GetGeometryRef()
249        resg=resg.GetBoundary()
250        tmp.SetGeometryDirectly(resg)
251        rgeometries+=[tmp]
252        geometry[i].Destroy()
253        i+=1
254    outputResult(conf,outputs["Result"],rgeometries)
255    return zoo.SERVICE_SUCCEEDED
256
257def CentroidPy(conf,inputs,outputs):
258    geometry=extractInputs(conf,inputs["InputPolygon"])
259    i=0
260    rgeometries=[]
261    while i < len(geometry):
262        tmp=geometry[i].Clone()
263        resg=geometry[i].GetGeometryRef()
264        if resg.GetGeometryType()!=3:
265            resg=resg.ConvexHull()
266        resg=resg.Centroid()
267        tmp.SetGeometryDirectly(resg)
268        rgeometries+=[tmp]
269        geometry[i].Destroy()
270        i+=1
271    outputResult(conf,outputs["Result"],rgeometries)
272    return zoo.SERVICE_SUCCEEDED
273
274def ConvexHullPy(conf,inputs,outputs):
275    geometry=extractInputs(conf,inputs["InputPolygon"])
276    i=0
277    rgeometries=[]
278    while i < len(geometry):
279        tmp=geometry[i].Clone()
280        resg=geometry[i].GetGeometryRef().ConvexHull()
281        tmp.SetGeometryDirectly(resg)
282        rgeometries+=[tmp]
283        geometry[i].Destroy()
284        i+=1
285    outputResult(conf,outputs["Result"],rgeometries)
286    return zoo.SERVICE_SUCCEEDED
287
288
289
290def EnvelopePy(conf,inputs,outputs):
291    print >> sys.stderr, inputs
292    try:
293        bdist=float(inputs["BufferDistance"]["value"])
294    except:
295        bdist=10
296    geometry=extractInputs(conf,inputs["InputPolygon"])
297    tmp=geometry[0].GetGeometryRef().GetEnvelope()
298    outputs["Result"]["value"]=str(tmp[0])+','+str(tmp[2])+','+str(tmp[1])+','+str(tmp[3])+','+'urn:ogc:def:crs:OGC:1.3:CRS84'
299    print >> sys.stderr,outputs["Result"]
300    return 3
301
302def UnionPy(conf,inputs,outputs):
303    geometry1=extractInputs(conf,inputs["InputEntity1"])
304    geometry2=extractInputs(conf,inputs["InputEntity2"])
305    rgeometries=[]
306    i=0
307    while i < len(geometry1):
308        j=0
309        while j < len(geometry2):
310            tmp=geometry1[i].Clone()
311            resg=geometry2[j].GetGeometryRef()
312            resg=resg.Union(geometry1[i].GetGeometryRef())
313            tmp.SetGeometryDirectly(resg)
314            if not(resg.IsEmpty()):
315                rgeometries+=[tmp]
316            j+=1
317        geometry1[i].Destroy()
318        i+=1
319    i=0
320    while i < len(geometry2):
321        geometry2[i].Destroy()
322        i+=1
323    outputResult(conf,outputs["Result"],rgeometries)
324    return 3
325
326def IntersectionPy(conf,inputs,outputs):
327
328    geometry1=extractInputs(conf,inputs["InputEntity1"])
329    geometry2=extractInputs(conf,inputs["InputEntity2"])
330
331    print >> sys.stderr,str(len(geometry1))+" "+str(len(geometry2))
332
333    rgeometries=[]
334    fids=[]
335    i=0
336    while i < len(geometry1):
337        j=0
338        while j < len(geometry2):
339            tmp=geometry2[j].Clone()
340            resg=geometry2[j].GetGeometryRef()
341            #resg=resg.Intersection(geometry1[i].GetGeometryRef())
342            resg=geometry1[i].GetGeometryRef().Intersection(resg)
343            tmp.SetGeometryDirectly(resg)
344            if resg is not None and not(resg.IsEmpty()) and fids.count(tmp.GetFID())==0:
345                rgeometries+=[tmp]
346                fids+=[tmp.GetFID()]
347            else:
348                tmp.Destroy()
349            j+=1
350        geometry1[i].Destroy()
351        i+=1
352    i=0
353    while i < len(geometry2):
354        geometry2[i].Destroy()
355        i+=1
356    outputResult(conf,outputs["Result"],rgeometries)
357    print >> sys.stderr,"/outputResult"
358    return 3
359
360def DifferencePy(conf,inputs,outputs):
361    geometry1=extractInputs(conf,inputs["InputEntity1"])
362    geometry2=extractInputs(conf,inputs["InputEntity2"])
363    rgeometries=[]
364    i=0
365    while i < len(geometry1):
366        j=0
367        while j < len(geometry2):
368            tmp=geometry2[j].Clone()
369            resg=geometry1[i].GetGeometryRef()
370            resg=resg.Difference(geometry2[i].GetGeometryRef())
371            tmp.SetGeometryDirectly(resg)
372            if not(resg.IsEmpty()):
373                rgeometries+=[tmp]
374            j+=1
375        geometry1[i].Destroy()
376        i+=1
377    i=0
378    while i < len(geometry2):
379        geometry2[i].Destroy()
380        i+=1
381    outputResult(conf,outputs["Result"],rgeometries)
382    return 3
383
384def SymDifferencePy(conf,inputs,outputs):
385    geometry1=extractInputs(conf,inputs["InputEntity1"])
386    geometry2=extractInputs(conf,inputs["InputEntity2"])
387    rgeometries=[]
388    i=0
389    while i < len(geometry1):
390        j=0
391        while j < len(geometry2):
392            tmp=geometry2[j].Clone()
393            resg=geometry1[i].GetGeometryRef()
394            resg=resg.SymmetricDifference(geometry2[i].GetGeometryRef())
395            tmp.SetGeometryDirectly(resg)
396            rgeometries+=[tmp]
397            j+=1
398        geometry1[i].Destroy()
399        i+=1
400    i=0
401    while i < len(geometry2):
402        geometry2[i].Destroy()
403        i+=1
404    outputResult(conf,outputs["Result"],rgeometries)
405    return 3
406
Note: See TracBrowser for help on using the repository browser.

Download now !

Search

Context Navigation

Main Navigation

ZOO Sponsors

http://www.zoo-project.org/trac/chrome/site/img/geolabs-logo.pnghttp://www.zoo-project.org/trac/chrome/site/img/neogeo-logo.png http://www.zoo-project.org/trac/chrome/site/img/apptech-logo.png http://www.zoo-project.org/trac/chrome/site/img/3liz-logo.png http://www.zoo-project.org/trac/chrome/site/img/gateway-logo.png

Become a sponsor !

Knowledge partners

http://www.zoo-project.org/trac/chrome/site/img/ocu-logo.png http://www.zoo-project.org/trac/chrome/site/img/gucas-logo.png http://www.zoo-project.org/trac/chrome/site/img/polimi-logo.png http://www.zoo-project.org/trac/chrome/site/img/fem-logo.png http://www.zoo-project.org/trac/chrome/site/img/supsi-logo.png http://www.zoo-project.org/trac/chrome/site/img/cumtb-logo.png

Become a knowledge partner

Related links

http://zoo-project.org/img/ogclogo.png http://zoo-project.org/img/osgeologo.png