From 0ab55c61f1daa7a8ea4f413211011599d28a8268 Mon Sep 17 00:00:00 2001 From: patricev Date: Tue, 5 Apr 2016 15:12:52 +0200 Subject: [PATCH] automatic recognition of water depth + bug fix --- dialogs/posttelemacpropertiesdialog.py | 14 +- libs/caduc/rasterinterpolator.py | 29 ++++ libs/post_telemac_pluginlayer.py | 60 +++++--- libs/posttelemac_util.py | 41 +++++- libs/posttelemac_util_flow.py | 16 +-- libs/posttelemac_util_graphtemp.py | 61 ++------ .../posttelemac_selafin_parser.py | 92 +++++++++++- toshape/posttelemac_util_extractshp.py | 84 ++++++++++- ui/properties.ui | 135 +++++++++++++++++- 9 files changed, 443 insertions(+), 89 deletions(-) create mode 100644 libs/caduc/rasterinterpolator.py diff --git a/dialogs/posttelemacpropertiesdialog.py b/dialogs/posttelemacpropertiesdialog.py index cc2f317..277fe89 100644 --- a/dialogs/posttelemacpropertiesdialog.py +++ b/dialogs/posttelemacpropertiesdialog.py @@ -225,6 +225,9 @@ def __init__(self, layer1, parent=None): self.pushButton_9.clicked.connect(self.set_utilcrs) self.pushButton_2.clicked.connect(self.create_shp_points) + #raster creation + self.pushButton_createraster.clicked.connect(self.postutils.rasterCreation) + #final action self.initTreewidgettoolsindextab() self.treeWidget_utils.expandAll() @@ -428,6 +431,7 @@ def open_def_variables(self, lst_param): self.populatecombobox_param() self.layer.updateSelafinValues() self.setTreeWidgetIndex(self.treeWidget_parameters,0,index) + def delete_def_variables(self): """ @@ -440,6 +444,11 @@ def delete_def_variables(self): self.layer.parametres[index:index+1] = [] self.populatecombobox_param() self.setTreeWidgetIndex(self.treeWidget_parameters,0,index-1) + #checkkeysparameter + self.layer.parametreh = None + self.layer.parametrevx = None + self.layer.parametrevy = None + self.layer.updateSelafinValues() #Display tools - contour - color ramp things *********************************************** @@ -699,9 +708,11 @@ def populatecombobox_param(self): Populate parameters comboboxes on dialog update """ self.comboBox_parametreschooser.clear() + self.comboBox_parametreschooser_2.clear() for i in range(len(self.layer.parametres)): temp1 = [str(self.layer.parametres[i][0])+" : "+str(self.layer.parametres[i][1])] self.comboBox_parametreschooser.addItems(temp1) + self.comboBox_parametreschooser_2.addItems(temp1) self.treeWidget_parameters.clear() itms = [] @@ -770,7 +781,8 @@ def initTreewidgettoolsindextab(self): [[-1,6],7, 'Max res' ] , [[7,0],8,'2shape contour' ], [[7,1],9,'2shape mesh' ], - [[7,2],10,'2shape point' ]] + [[7,2],10,'2shape point' ], + [[8,0],11,'Raster creation']] def changepannelutils(self): diff --git a/libs/caduc/rasterinterpolator.py b/libs/caduc/rasterinterpolator.py new file mode 100644 index 0000000..bd81bc0 --- /dev/null +++ b/libs/caduc/rasterinterpolator.py @@ -0,0 +1,29 @@ +""" +ld=qgis.analysis.QgsIDWInterpolator.LayerData +ld.vectorLayer=iface.activeLayer() + +ld.zCoordInterpolation=0 +ld.interpolationAttribute=5 +ld.mInputType=0 + +#print str( qgis.analysis.QgsIDWInterpolator.layerData() ) +""" + +ld1 = qgis.analysis.QgsInterpolator.LayerData() +ld1.vectorLayer=iface.activeLayer() +ld1.zCoordInterpolation=0 +ld1.interpolationAttribute=5 +ld1.mInputType=0 + +#print dir(qgis.analysis.QgsIDWInterpolator) +#print dir(qgis.analysis.QgsInterpolator.LayerData) + +itp=qgis.analysis.QgsIDWInterpolator([ld1]) + + +rect = iface.mapCanvas().extent() +ncol = 10 +res = 1.0 +test = qgis.analysis.QgsGridFileWriter(itp,'c:/test1.asc',rect,ncol,ncol,res,res) +test.writeFile(True) + diff --git a/libs/post_telemac_pluginlayer.py b/libs/post_telemac_pluginlayer.py index 4ec9481..ea10874 100644 --- a/libs/post_telemac_pluginlayer.py +++ b/libs/post_telemac_pluginlayer.py @@ -32,7 +32,7 @@ #import scipy #from scipy.spatial import cKDTree #import matplotlib -#from matplotlib import tri +from matplotlib import tri from matplotlib import colors #import PyQT from PyQt4.QtCore import * @@ -251,24 +251,43 @@ def initSelafinParameters(self): if len(self.parametrestoload)>0: #case of virtual parameters when loadin a selafin layer for param in self.parametrestoload: self.parametres.append([len(self.parametres),param[1],param[2],len(self.parametres)]) - try: #load velocity parameters - self.parametrevx = self.propertiesdialog.postutils.getParameterName("VITESSEU")[0] - self.parametrevy = self.propertiesdialog.postutils.getParameterName("VITESSEV")[0] - self.propertiesdialog.tab_velocity.setEnabled(True) - for widget in self.propertiesdialog.tab_velocity.children(): - widget.setEnabled(True) - for widget in self.propertiesdialog.groupBox_schowvel.children(): - widget.setEnabled(True) - self.propertiesdialog.groupBox_schowvel.setChecked(True) - self.propertiesdialog.groupBox_schowvel.setChecked(False) - except Exception, e: - self.propertiesdialog.tab_velocity.setEnabled(False) - #TODO : disable utils dependant on velocity (flow, show velocity) - try: #load water depth parameters - self.parametreh = self.propertiesdialog.postutils.getParameterName("HAUTEUR")[0] - except Exception, e: - pass - #TODO : disable utils dependant on velocity (flow, show velocity) + + self.identifyKeysParameters() + + + def identifyKeysParameters(self): + #load velocity parameters + if (self.parametrevx == None and self.parametrevy == None) : + + if self.propertiesdialog.postutils.getParameterName("VITESSEU") == None and self.propertiesdialog.postutils.getParameterName("VITESSEV") == None: + self.propertiesdialog.tab_velocity.setEnabled(False) + else: + self.parametrevx = self.propertiesdialog.postutils.getParameterName("VITESSEU")[0] + self.parametrevy = self.propertiesdialog.postutils.getParameterName("VITESSEV")[0] + self.propertiesdialog.tab_velocity.setEnabled(True) + for widget in self.propertiesdialog.tab_velocity.children(): + widget.setEnabled(True) + for widget in self.propertiesdialog.groupBox_schowvel.children(): + widget.setEnabled(True) + self.propertiesdialog.groupBox_schowvel.setChecked(True) + self.propertiesdialog.groupBox_schowvel.setChecked(False) + + #load water depth parameters + if self.parametreh == None : + if self.propertiesdialog.postutils.getParameterName("HAUTEUR") == None: + paramfreesurface = self.propertiesdialog.postutils.getParameterName("SURFACELIBRE")[0] + parambottom = self.propertiesdialog.postutils.getParameterName("BATHYMETRIE")[0] + self.parametreh = len(self.parametres) + self.parametres.append([len(self.parametres),"HAUTEUR D'EAU",'V'+str(paramfreesurface)+' - V'+str(parambottom),len(self.parametres)]) + else: + self.parametreh = self.propertiesdialog.postutils.getParameterName("HAUTEUR")[0] + + if self.parametreh != None and self.parametrevx != None and self.parametrevy != None: + self.propertiesdialog.page_flow.setEnabled(True) + else: + self.propertiesdialog.page_flow.setEnabled(False) + + def clearParameters(self): self.param_displayed = None @@ -304,6 +323,7 @@ def updateSelafinValues(self, onlyparamtimeunchanged = -1 ): It emits a signal, because this method is redirected to a different method in comparetool when tool "compare" is activated """ + self.identifyKeysParameters() self.updatevalue.emit(onlyparamtimeunchanged) updatevalue = pyqtSignal(int) @@ -576,7 +596,7 @@ def identify(self,qgspointfromcanvas,multiparam = False): return (True,d) def initTriinterpolator(self): - self.triinterp = [tri.LinearTriInterpolator(self.hydrauparser.triangulation, self.values[i]) for i in range(len(self.parametres))] + self.triinterp = [matplotlib.tri.LinearTriInterpolator(self.hydrauparser.triangulation, self.values[i]) for i in range(len(self.parametres))] #**************************************************************************************************** #************Method for saving/loading project with selafinlayer file*********************************** diff --git a/libs/posttelemac_util.py b/libs/posttelemac_util.py index 29872af..a0a593d 100644 --- a/libs/posttelemac_util.py +++ b/libs/posttelemac_util.py @@ -27,6 +27,7 @@ from PyQt4.QtCore import * from PyQt4.QtGui import * from PyQt4.QtCore import SIGNAL, Qt +import qgis #import numpy import numpy as np #import time @@ -342,7 +343,7 @@ def launchThread(self,geom): def workerFinished(self,list1,list2,list3 = None): - + if self.graphtodo ==0: ax = self.selafinlayer.propertiesdialog.ax if not self.selafinlayer.propertiesdialog.checkBox.isChecked(): @@ -634,7 +635,43 @@ def workerFinishedHillshade(self,strpath): #self.selafinlayer.propertiesdialog.textBrowser_main.append(ctime() + " - "+ str(os.path.basename(strpath).split('.')[0]) + self.tr(" created")) self.selafinlayer.propertiesdialog.normalMessage(str(os.path.basename(strpath).split('.')[0]) + self.tr(" created")) - + def rasterCreation(self): + #print 'ok' + #Points + layerstring = self.selafinlayer.hydraufilepath+'[' + str(self.selafinlayer.time_displayed) + ']' + print layerstring + layerpoint = QgsVectorLayer(layerstring, 'test', "ogr") + #Layerpoint for interpolation + ld1 = qgis.analysis.QgsInterpolator.LayerData() + ld1.vectorLayer=layerpoint + ld1.zCoordInterpolation=False + ld1.interpolationAttribute=self.selafinlayer.propertiesdialog.comboBox_parametreschooser_2.currentIndex() + ld1.mInputType=0 + + #interpolator + #itp=qgis.analysis.QgsIDWInterpolator([ld1]) + itp=qgis.analysis.QgsTINInterpolator([ld1]) + + if self.selafinlayer.propertiesdialog.comboBox_rasterextent.currentIndex() == 0 : + rect = iface.mapCanvas().extent() + elif self.selafinlayer.propertiesdialog.comboBox_rasterextent.currentIndex() == 1 : + rect = layerpoint.extent() + + #raster creation + res = self.selafinlayer.propertiesdialog.spinBox_rastercellsize.value() + ncol = int( ( rect.xMaximum() - rect.xMinimum() ) / res ) + print str(ncol) + ncol = 300 + rasterfilepath = os.path.join(os.path.dirname(self.selafinlayer.hydraufilepath),'raster2.asc') + print rasterfilepath + test = qgis.analysis.QgsGridFileWriter(itp,rasterfilepath,rect,ncol,ncol,res,res) + print 'ok1' + test.writeFile(False) + + + + + #**************************************************************************************************** diff --git a/libs/posttelemac_util_flow.py b/libs/posttelemac_util_flow.py index 6bb4aa7..7ec9b39 100644 --- a/libs/posttelemac_util_flow.py +++ b/libs/posttelemac_util_flow.py @@ -94,12 +94,12 @@ def computeFlowMain(self): try: if i==0: #init try: - h2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameterh])[0][0]) + h2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameterh],self.selafinlayer.parametres)[0][0]) except Exception , e : self.status.emit('method 011 : ' + str(e)) - uv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameteruv])[0][0]) + uv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameteruv],self.selafinlayer.parametres)[0][0]) uv2 = np.array([[value,0.0] for value in uv2]) - vv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parametervv])[0][0]) + vv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parametervv],self.selafinlayer.parametres)[0][0]) vv2 = np.array([[0.0,value] for value in vv2]) v2vect = uv2 + vv2 #xy2 = [self.slf.MESHX[elem],self.slf.MESHY[elem]] @@ -108,10 +108,10 @@ def computeFlowMain(self): h1 = h2 v1vect = v2vect xy1 = xy2 - h2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameterh])[0][0]) - uv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameteruv])[0][0]) + h2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameterh],self.selafinlayer.parametres)[0][0]) + uv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parameteruv],self.selafinlayer.parametres)[0][0]) uv2 = np.array([[value,0.0] for value in uv2]) - vv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parametervv])[0][0]) + vv2 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([elem + 1],[parametervv],self.selafinlayer.parametres)[0][0]) vv2 = np.array([[0.0,value] for value in vv2]) v2vect = uv2 + vv2 #xy2 = [self.slf.MESHX[elem],self.slf.MESHY[elem]] @@ -384,8 +384,8 @@ def computeFlowBetweenPoints(self,xy1,h1,v1vect,xy2,h2,v2vect): def valuebetweenEdges(self,xy,edges,param): xytemp = np.array(xy) - h11 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([edges[0] + 1],[param])[0][0]) #getseries begins at 1 - h12 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([edges[1] + 1 ],[param])[0][0]) + h11 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([edges[0] + 1],[param],self.selafinlayer.parametres)[0][0]) #getseries begins at 1 + h12 = np.array(self.selafinlayer.hydrauparser.getTimeSerie([edges[1] + 1 ],[param],self.selafinlayer.parametres)[0][0]) """ e1 = np.array([self.selafinlayer.slf.MESHX[edges[0]],self.selafinlayer.slf.MESHY[edges[0]]]) e2 = np.array([self.selafinlayer.slf.MESHX[edges[1]],self.selafinlayer.slf.MESHY[edges[1]]]) diff --git a/libs/posttelemac_util_graphtemp.py b/libs/posttelemac_util_graphtemp.py index 3ca22f9..5d6e547 100644 --- a/libs/posttelemac_util_graphtemp.py +++ b/libs/posttelemac_util_graphtemp.py @@ -57,14 +57,21 @@ def createGraphTemp(self): abscisse = self.selafinlayer.hydrauparser.getTimes().tolist() param=self.selafinlayer.propertiesdialog.comboBox_parametreschooser.currentIndex() - #param = self.selafinlayer.propertiesdialog.getTreeWidgetSelectedIndex(self.selafinlayer.propertiesdialog.treeWidget_parameters)[1] - if self.selafinlayer.parametres[param][2]: - dico = self.getDico(self.selafinlayer.parametres[param][2], self.selafinlayer.parametres, self.selafinlayer.values,enumpoint) - tempordonees = eval(self.selafinlayer.parametres[param][2],{}, dico) + + if self.compare : + triangles,numpointsfinal,pointsfinal,coef = self.selafinlayer.propertiesdialog.postutils.compareprocess.hydrauparsercompared.getInterpFactorInTriangleFromPoint([x],[y]) + self.status.emit(str(triangles)+' ' +str(numpointsfinal)+' ' +str(pointsfinal)+' ' +str(coef)) + layer2serie = 0 + print str(numpointsfinal[0]) + for i, numpoint in enumerate(numpointsfinal[0]): + #layer2serie += float(coef[0][i]) * self.selafinlayer.propertiesdialog.postutils.compareprocess.hydrauparsercompared.getTimeSerie([numpoint],[self.selafinlayer.parametres[param[0]][3]],self.selafinlayer.parametres) + layer2serie += float(coef[0][i]) * self.selafinlayer.propertiesdialog.postutils.compareprocess.hydrauparsercompared.getTimeSerie([numpoint +1],[self.selafinlayer.parametres[param][3]],self.selafinlayer.parametres) + print 'ok1' + layer1serie = self.selafinlayer.hydrauparser.getTimeSerie([enumpoint + 1],[param],self.selafinlayer.parametres) + tempordonees = layer2serie - layer1serie else: - #tempordonees = self.selafinlayer.slf.getSERIES([enumpoint + 1],[param],False) #points in getseries begin with 1 - #tempordonees = self.selafinlayer.hydrauparser.getTimeSerie([enumpoint + 1],[param]) #points in getseries begin with 1 - tempordonees = self.getGraphTempSeries([enumpoint + 1],[param]) #points in getseries bein with 1 + tempordonees = self.selafinlayer.hydrauparser.getTimeSerie([enumpoint + 1],[param],self.selafinlayer.parametres) + ordonnees = tempordonees[0][0].tolist() list1.append(abscisse) list2.append(ordonnees) @@ -72,47 +79,7 @@ def createGraphTemp(self): except Exception, e: self.status.emit(str(e)) self.finished.emit([],[]) - - - def getDico(self,expr, parametres, values,enumpoint): - dico = {} - try: - - dico['sin'] = sin - dico['cos'] = cos - dico['abs'] = abs - dico['int'] = int - dico['if_then_else'] = self.selafinlayer.if_then_else - a = 'V{}' - nb_var = len(values) - i = 0 - num_var = 0 - while num_var < nb_var: - if not parametres[i][2]: - #dico[a.format(i)] = self.selafinlayer.hydrauparser.getTimeSerie([enumpoint + 1],[i]) - dico[a.format(i)] = self.getGraphTempSeries([enumpoint + 1],[i]) - num_var += 1 - i += 1 - except Exception, e: - print str(e) - return dico - - def getGraphTempSeries(self,num,param): - if self.compare : - x,y = self.selafinlayer.hydrauparser.getXYFromNumPoint(num)[0] - triangles,numpointsfinal,pointsfinal,coef = self.selafinlayer.propertiesdialog.postutils.compareprocess.hydrauparsercompared.getInterpFactorInTriangleFromPoint([x],[y]) - layer2serie = 0 - for i, numpoint in enumerate(numpointsfinal[0]): - layer2serie += float(coef[0][i]) * self.selafinlayer.propertiesdialog.postutils.compareprocess.hydrauparsercompared.getTimeSerie([numpoint],[self.selafinlayer.parametres[param[0]][3]]) - layer1serie = self.selafinlayer.hydrauparser.getTimeSerie(num,param) - return layer2serie - layer1serie - else: - return self.selafinlayer.hydrauparser.getTimeSerie(num,param) - - - - progress = QtCore.pyqtSignal(int) status = QtCore.pyqtSignal(str) error = QtCore.pyqtSignal(str) diff --git a/posttelemacparsers/posttelemac_selafin_parser.py b/posttelemacparsers/posttelemac_selafin_parser.py index 29d2b63..b59f95c 100644 --- a/posttelemacparsers/posttelemac_selafin_parser.py +++ b/posttelemacparsers/posttelemac_selafin_parser.py @@ -22,6 +22,7 @@ """ import matplotlib from matplotlib import tri +from numpy import * import numpy as np from scipy.spatial import cKDTree from ..libs_telemac.parsers.parserSELAFIN import SELAFIN @@ -60,8 +61,27 @@ def getValues(self,time): return self.hydraufile.getVALUES(time) - def getTimeSerie(self,arraynumpoint,arrayparam): - return self.hydraufile.getSERIES(arraynumpoint,arrayparam,False) + def getTimeSerie(self,arraynumpoint,arrayparam,layerparametres = None): + """ + Warning : point index begin at 1 + """ + if False: + return self.hydraufile.getSERIES(arraynumpoint,arrayparam,False) + else: + #print 'arraypoint ' + str(arraynumpoint) + result = [] + try: + for param in arrayparam: + if layerparametres != None and layerparametres[param][2]: + dico = self.getDico(layerparametres[param][2], layerparametres,arraynumpoint) + tempordonees = eval(layerparametres[param][2],{}, dico) + result.append(tempordonees[0]) + else: + tempordonees = self.hydraufile.getSERIES(arraynumpoint,[param],False) + result.append(tempordonees[0]) + except Exception, e: + print 'getserie ' + str(e) + return np.array(result) def getMesh(self): return (self.hydraufile.MESHX, self.hydraufile.MESHY) @@ -76,7 +96,75 @@ def getIkle(self): def getTimes(self): return self.hydraufile.tags["times"] + + #Others method - don't touch it + + def getDico(self,expr, parametres,enumpoint): + dico = {} + try: + dico['sin'] = sin + dico['cos'] = cos + dico['abs'] = abs + dico['int'] = int + dico['if_then_else'] = self.if_then_else + a = 'V{}' + #nb_var = len(values) + nb_var = len( self.getValues(0) ) + i = 0 + num_var = 0 + while num_var < nb_var: + if not parametres[i][2]: + dico[a.format(i)] = self.getTimeSerie(enumpoint,[i]) + num_var += 1 + i += 1 + except Exception, e: + print 'getdico ' + str(e) + return dico + + def if_then_else(self,ifstat,true1,false1): + """ + Used for calculation of virtual parameters + """ + #condition + if isinstance(ifstat,np.ndarray): + var2 = np.zeros(ifstat.shape) + temp1 = np.where(ifstat) + elif isinstance(ifstat,str): + val = eval(ifstat,{"__builtins__":None}, self.dico) + var2 = np.zeros(val.shape) + temp1 = np.where(val) + #True + if isinstance(true1,np.ndarray): + var2[temp1] = true1[temp1] + elif isinstance(true1, numbers.Number): + var2[temp1] = float(true1) + else: + pass + #False + mask = np.ones(var2.shape, np.bool) + mask[temp1] = 0 + if isinstance(false1,np.ndarray): + var2[mask] = false1[mask] + elif isinstance(false1, numbers.Number): + var2[mask] = float(false1) + else: + pass + return var2 + """ + def getGraphTempSeries(self,num,param): + if self.compare : + x,y = self.getXYFromNumPoint(num)[0] + triangles,numpointsfinal,pointsfinal,coef = self.getInterpFactorInTriangleFromPoint([x],[y]) + layer2serie = 0 + for i, numpoint in enumerate(numpointsfinal[0]): + layer2serie += float(coef[0][i]) * self.selafinlayer.propertiesdialog.postutils.compareprocess.hydrauparsercompared.getTimeSerie([numpoint],[self.selafinlayer.parametres[param[0]][3]]) + layer1serie = self.selafinlayer.hydrauparser.getTimeSerie(num,param) + return layer2serie - layer1serie + else: + return self.selafinlayer.hydrauparser.getTimeSerie(num,param) + """ + def getXYFromNumPoint(self,arraynumpoint): meshx,meshy = self.getMesh() diff --git a/toshape/posttelemac_util_extractshp.py b/toshape/posttelemac_util_extractshp.py index b701d0e..d5dd65f 100644 --- a/toshape/posttelemac_util_extractshp.py +++ b/toshape/posttelemac_util_extractshp.py @@ -23,7 +23,10 @@ from time import ctime import math from os import path -#from shapely.geometry import Polygon +#shapely +from shapely import * +from shapely.geometry import Polygon +from shapely.wkb import loads import sys import os.path """ @@ -353,11 +356,23 @@ def get_outerinner(self,geom1): def InsertRinginFeature(self,f1,allfeatures2,vlOuterTempIndex,lvltemp1,counttotal): try: #Correction des erreurs de geometrie des outers - if len(f1.geometry().validateGeometry())!=0: - f1geom = f1.geometry().buffer(0.01,5) - if f1geom.area() < f1.geometry().area(): - f1geom = f1.geometry() - self.writeOutput(ctime() + ' - Warning : geometry '+str(f1.id())+' not valid before inserting rings') + if len(f1.geometry().validateGeometry()) != 0: + if True: + f1geom = f1.geometry().buffer(0.01,5) + if f1geom.area() < f1.geometry().area(): + f1geom = f1.geometry() + self.writeOutput(ctime() + ' - Warning : geometry '+str(f1.id())+' not valid before inserting rings') + else: + geomshpapely = loads(f1.geometry().asWkb()) + resulttemp = self.repairPolygon(geomshpapely) + if resulttemp != None: + geom3 = [QgsPoint(point[0],point[1]) for point in list( resulttemp.exterior.coords )] + geom2 = QgsGeometry.fromPolygon([geom3]) + f1geom = geom2 + self.writeOutput(ctime() + ' - geometry correction') + else: + f1geom = f1.geometry() + else: f1geom = f1.geometry() @@ -372,7 +387,8 @@ def InsertRinginFeature(self,f1,allfeatures2,vlOuterTempIndex,lvltemp1,counttota for id in ids: f2geom = allfeatures2[id].geometry() if len(f2geom.validateGeometry())!=0: - f2geom= f2geom.buffer(0.01,5) + f2geom= f2geom.buffer(0.00,5) + tab.append([f2geom.area(),f2geom]) if len(tab)>0: tab.sort(reverse = True) @@ -420,6 +436,60 @@ def do_ring(self,geom3): ring.append(QgsPoint(polygon[i][0],polygon[i][1])) ring.append(QgsPoint(polygon[0][0],polygon[0][1])) return ring + + + def repairPolygon(self,geometry): + buffer_worker = True + try: + geometry = geometry.buffer(0) + except : + buffer_worker = False + + if buffer_worker: + return geometry + + polygons = [] + if geometry.geom_type == "Polygon": + polygons.append(geometry) + elif geometry.geom_type == "MultiPolygon": + polygons.extend(geometry.geoms) + + + fixed_polygons = [] + for n, polygon in enumerate(polygons): + if not self.linear_ring_is_valid(polygon.exterior): + continue #"unable to fix" + + interiors = [] + for ring in polygon.interiors: + if self.linear_ring_is_valid(ring): + interiors.append(ring) + + fixed_polygon = shapely.geometry.Polygon(polygon.exterior, interiors) + + try: + fixed_polygon = fixed_polygon.buffer(0) + except : + continue + + if fixed_polygon.geom_type == "Polygon": + fixed_polygons.append(fixed_polygon) + elif fixed_polygon.geom_type == "MultiPolygon": + fixed_polygons.extend(fixed_polygon.geoms) + + if len(fixed_polygons) > 0 : + return shapely.geometry.MultiPolygon(fixed_polygons) + else: + return None + + def linear_ring_is_valid(self,ring): + points = set() + for x,y in ring.coords: + points.add((x,y)) + if len(points) < 3: + return False + else: + return True progress = QtCore.pyqtSignal(int) diff --git a/ui/properties.ui b/ui/properties.ui index 0e8eea3..a03fd69 100644 --- a/ui/properties.ui +++ b/ui/properties.ui @@ -7,7 +7,7 @@ 0 0 554 - 641 + 711 @@ -1491,6 +1491,16 @@ max .res + + + Rasters + + + + Raster creation + + + @@ -1841,7 +1851,10 @@ max .res - + + + false + @@ -2751,6 +2764,124 @@ max .res + + + + + + + 0 + 0 + + + + + 16777215 + 90 + + + + + + + + 20 + + + + + + 0 + 0 + + + + + + + + + 0 + 0 + + + + Parameter + + + + + + + + + + GroupBox + + + + + + Cell size (0 : automatic) + + + + + + + 999999 + + + 2 + + + + + + + Extent + + + + + + + + Current view + + + + + Full Extent + + + + + + + + + + + Create raster + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + +