From 23523aa05de648795e084017e35664720236fc19 Mon Sep 17 00:00:00 2001 From: lennepkade Date: Sun, 12 Jun 2016 10:42:58 -0300 Subject: [PATCH] v1.1 --- dzetsaka.py | 53 ++++++++------- metadata.txt | 2 +- scripts/filters.py | 120 ++++++++++++++++++--------------- scripts/function_dataraster.py | 2 +- scripts/mainfunction.py | 5 +- ui/dzetsaka_dock.py | 2 +- 6 files changed, 100 insertions(+), 84 deletions(-) diff --git a/dzetsaka.py b/dzetsaka.py index b352a88..909f5d4 100644 --- a/dzetsaka.py +++ b/dzetsaka.py @@ -122,7 +122,6 @@ def __init__(self, iface): self.dockwidget.outMatrix.clear() self.dockwidget.checkOutMatrix.clicked.connect(self.checkbox_state) - ## init fields self.dockwidget.inField.clear() self.dockwidget.inShape.currentIndexChanged[int].connect(self.onChangedLayer) @@ -134,11 +133,14 @@ def __init__(self, iface): fields = provider.fields() listFieldNames = [field.name() for field in fields] self.dockwidget.inField.addItems(listFieldNames) - + + ## + + self.dockwidget.setMaximumHeight(400) ## let's run the classification ! self.dockwidget.performMagic.clicked.connect(self.runMagic) - + def onChangedLayer(self): """!@brief If active layer is changed, change column combobox""" # We clear combobox @@ -156,8 +158,11 @@ def onChangedLayer(self): def select_output_file(self): """!@brief Select file to save, and gives the right extension if the user don't put it""" sender = self.sender() - - fileName = QFileDialog.getSaveFileName(self.dockwidget, "Select output file","","TIF (*.tif)") + + if sender == self.historicalmap.outShpButton: + fileName = QFileDialog.getSaveFileName(self.dockwidget, "Select output file","","SHP (*.shp)") + else: + fileName = QFileDialog.getSaveFileName(self.dockwidget, "Select output file","","TIF (*.tif)") if not fileName: return @@ -180,13 +185,14 @@ def select_output_file(self): self.historicalmap.outShp.setText(fileName+'.shp') else: self.historicalmap.outShp.setText(fileName+fileExtension) - - if sender == self.filters_dock.outRasterButton: - if fileExtension!='.tif': - self.filters_dock.outRaster.setText(fileName+'.tif') + try: + if sender == self.filters_dock.outRasterButton: + if fileExtension!='.tif': + self.filters_dock.outRaster.setText(fileName+'.tif') else: self.filters_dock.outRaster.setText(fileName+fileExtension) - + except: + print('ko') def checkbox_state(self): sender=self.sender() @@ -510,7 +516,7 @@ def loadFilters(self): self.filters_dock.inFilter.addItems(filtersList) # make choice - if self.sender() == self.menu.filterErosion: + if self.sender() == self.menu.filterOpening: self.filters_dock.inFilter.setCurrentIndex(0) elif self.sender() == self.menu.filterClosing: self.filters_dock.inFilter.setCurrentIndex(1) @@ -559,7 +565,7 @@ def runFilter(self): inFilterIter = self.filters_dock.inFilterIter.value() if self.filters_dock.outRaster.text()=='': - self.outRaster= tempfile.mkstemp('.tif')[1] + self.outRaster= tempfile.mktemp('.tif') else: self.outRaster= self.dockwidget.outRaster.text() @@ -605,7 +611,7 @@ def runHistoricalMapStep1(self): inRaster=inRaster.dataProvider().dataSourceUri() if self.historicalmap.outRaster.text()=='': - outRaster = tempfile.mkstemp('.tif')[1] + outRaster = tempfile.mktemp('.tif') else: outRaster = self.historicalmap.outRaster.text() @@ -633,20 +639,21 @@ def runHistoricalMapStep2(self): inRaster=inRaster.dataProvider().dataSourceUri() if self.historicalmap.outShp.text()=='': - outShp = tempfile.mkstemp('.shp')[1] + outShp = tempfile.mktemp('.shp') else: outShp = self.historicalmap.outShp.text() sieveSize = self.historicalmap.sieveSize.value() inClassNumber = self.historicalmap.classNumber.value() - try: - worker = filtersFunction() - outRaster = worker.historicalMapPostRaster(inRaster,sieveSize,inClassNumber,outShp) - outShp = worker.historicalMapPostVector(outRaster,outShp) - self.iface.addVectorLayer(outShp,os.path.splitext(os.path.basename(outShp))[0],'ogr') - except: - QtGui.QMessageBox.warning(self, 'Problem with Processing', 'Did you activate \'Processing\' plugin in Qgis ?', QtGui.QMessageBox.Ok) + + worker = filtersFunction() + outShp = worker.historicalMapPostClass(inRaster,sieveSize,inClassNumber,outShp) + + #outShp = worker.historicalMapPostVector(outRaster,outShp) + self.iface.addVectorLayer(outShp,os.path.splitext(os.path.basename(outShp))[0],'ogr') + #except: + # QtGui.QMessageBox.warning(self, 'Problem with Processing', 'Did you activate \'Processing\' plugin in Qgis ?', QtGui.QMessageBox.Ok) @@ -678,7 +685,7 @@ def runMagic(self): else: # create temp if not output raster if self.dockwidget.outRaster.text()=='': - outRaster= tempfile.mkstemp()[1] + outRaster= tempfile.mktemp('.tif') else: outRaster= self.dockwidget.outRaster.text() @@ -707,7 +714,7 @@ def runMagic(self): # Perform training else: if self.dockwidget.outModel.text()=='': - model=tempfile.mkstemp()[1] + model=tempfile.mktemp('.'+str(model)) else: model=self.dockwidget.outModel.text() diff --git a/metadata.txt b/metadata.txt index 4e7edef..c8ab973 100644 --- a/metadata.txt +++ b/metadata.txt @@ -23,7 +23,7 @@ repository=http://www.github.com/lennepkade/dzetsaka # Recommended items: # Uncomment the following line and add your changelog: -changelog= filters for image (closing, median...), and historical map process (for automatic vectorizing forest from old maps) +changelog= filters for image (closing, median...), and new historical map process (for automatic vectorizing forest from old maps) # Tags are comma separated with spaces allowed tags=classification,automatic,gaussian,mixture,model diff --git a/scripts/filters.py b/scripts/filters.py index e56ccad..9785996 100644 --- a/scripts/filters.py +++ b/scripts/filters.py @@ -100,47 +100,49 @@ def historicalMapFilter(self,inImage,outRaster,inShapeGrey,inShapeMedian,iterMed except: filterProgress.reset() - def historicalMapPostRaster(self,inRaster,sieveSize,inClassNumber,outShp): - import processing - - rasterTemp = tempfile.mkstemp('.tif')[1] - - processing.runalg('gdalogr:sieve',inRaster,sieveSize,0,rasterTemp) + def historicalMapPostClass(self,inRaster,sieveSize,inClassNumber,outShp): - ### remove unwanted classe + Progress=progressBar(' Post-classification...',10) + + rasterTemp = tempfile.mktemp('.tif') + + datasrc = gdal.Open(inRaster) + srcband = datasrc.GetRasterBand(1) + data,im=dataraster.open_data_band(inRaster) - data,im=dataraster.open_data_band(rasterTemp) + drv = gdal.GetDriverByName('GTiff') + dst_ds = drv.Create(rasterTemp,datasrc.RasterXSize,datasrc.RasterXSize,1,srcband.DataType) + dst_ds.SetGeoTransform(datasrc.GetGeoTransform()) + dst_ds.SetProjection(datasrc.GetProjection()) + + dstband=dst_ds.GetRasterBand(1) - # get proj,geo and dimension (d) from data - proj = data.GetProjection() - geo = data.GetGeoTransform() + def sieve(srcband,dstband,sieveSize): + gdal.SieveFilter(srcband,None,dstband,sieveSize) - outFile=dataraster.create_empty_tiff(rasterTemp,im,1,geo,proj) - - try: - temp = data.GetRasterBand(1).ReadAsArray() + sieve(srcband,dstband,sieveSize) + Progress.addStep(2) + dst_ds =None + + def reclass(rasterTemp,inClassNumber): + dst_ds = gdal.Open(rasterTemp,gdal.GA_Update) + srcband = dst_ds.GetRasterBand(1).ReadAsArray() # All data which is not forest is set to 0, so we fill all for the forest only, because it's a binary fill holes. # Set selected class as 1 - temp[temp!=inClassNumber]=0 - temp[temp==inClassNumber]=1 - #temp = ndimage.median_filter(temp,size=(3,3)).astype(int) - except: - QgsMessageLog.logMessage("Cannot sieve") - - out=outFile.GetRasterBand(1) - out.WriteArray(temp) - out.FlushCache() - temp = None - # Cleaning outFile or vectorizing doesn't work - outFile= None - - return rasterTemp + srcband[srcband !=inClassNumber]=0 + srcband[srcband ==inClassNumber]=1 + #QgsMessageLog.logMessage("Cannot sieve") + + dst_ds.GetRasterBand(1).WriteArray(srcband) + dst_ds.FlushCache() + return rasterTemp - def historicalMapPostVector(self,inRaster,outShp): + rasterTemp = reclass(rasterTemp,inClassNumber) + Progress.addStep(1) - historicalProgress=progressBar('Vectorizing...',2) - try: - sourceRaster = gdal.Open(inRaster) + def polygonize(rasterTemp,outShp): + + sourceRaster = gdal.Open(rasterTemp) band = sourceRaster.GetRasterBand(1) driver = ogr.GetDriverByName("ESRI Shapefile") # If shapefile already exist, delete it @@ -157,33 +159,38 @@ def historicalMapPostVector(self,inRaster,outShp): newField = ogr.FieldDefn('Class', ogr.OFTInteger) outLayer.CreateField(newField) + print('vectorizing') gdal.Polygonize(band, None,outLayer, 0,[],callback=None) + print('vectorized') outDatasource.Destroy() sourceRaster=None + band=None + # QgsMessageLog.logMessage("Cannot vectorize "+rasterTemp) + + ioShpFile = ogr.Open(outShp, update = 1) - except: - QgsMessageLog.logMessage("Cannot vectorize "+rasterTemp) - - ioShpFile = ogr.Open(outShp, update = 1) - - historicalProgress.addStep() - - lyr = ioShpFile.GetLayerByIndex(0) - lyr.ResetReading() - - for i in lyr: - # feat = lyr.GetFeature(i) - lyr.SetFeature(i) - lyr.SetFeature(i) - # if area is less than inMinSize or if it isn't forest, remove polygon - if i.GetField('Class')!=1: - lyr.DeleteFeature(i.GetFID()) - ioShpFile.Destroy() + lyr = ioShpFile.GetLayerByIndex(0) + nbFeatures = lyr.GetFeatureCount() + lyr.ResetReading() + + for i in lyr: + lyr.SetFeature(i) + # if area is less than inMinSize or if it isn't forest, remove polygon + if i.GetField('Class')!=1: + lyr.DeleteFeature(i.GetFID()) + ioShpFile.Destroy() + + #historicalProgress.reset() + + return outShp - historicalProgress.reset() + outShp = polygonize(rasterTemp,outShp) + Progress.addStep(7) + Progress.reset() return outShp - + + def filters(self,inImage,outRaster,inFilter,inFilterSize,inFilterIter): # open data with Gdal self.processed = 0 @@ -289,13 +296,14 @@ def kill(self): if __name__=="__main__": - inRaster = '/home/lennepkade/Bureau/datapag/02-Results/02-Data/spot/pansharp-Spot7_arvi.tif' + inRaster = '/home/lennepkade/Bureau/datapag/02-Results/02-Data/spot/pred.tif' outRaster = '/home/lennepkade/Bureau/datapag/02-Results/02-Data/spot/closing.tif' inFilter = 'Gaussian' inFilterSize = 10 inFilterIter = 1 worker=filtersFunction() - worker.filters(inRaster,outRaster,inFilter,10,1) - + #worker.filters(inRaster,outRaster,inFilter,10,1) + worker.historicalMapPostClass(inRaster,40,2,'/home/lennepkade/Bureau/datapag/02-Results/02-Data/spot/closing.shp') + #worker.historicalMapReclass('/tmp/nana2.tif',2) print(outRaster) diff --git a/scripts/function_dataraster.py b/scripts/function_dataraster.py index 6ce0b1f..915debf 100644 --- a/scripts/function_dataraster.py +++ b/scripts/function_dataraster.py @@ -15,7 +15,7 @@ def open_data_band(filename): im : empty table with right dimension (array) """ - data = gdal.Open(filename,gdal.GA_ReadOnly) + data = gdal.Open(filename,gdal.GA_Update) if data is None: print 'Impossible to open '+filename exit() diff --git a/scripts/mainfunction.py b/scripts/mainfunction.py index 0da9e5d..75f809f 100644 --- a/scripts/mainfunction.py +++ b/scripts/mainfunction.py @@ -458,11 +458,12 @@ def __init__(self,inMsg=' Loading...',inMaxStep=1): # set Maximum for progressBar prgBar.setMaximum(inMaxStep) - def addStep(self): + def addStep(self,step=1): """!@brief Add a step to the progressBar addStep() simply add +1 to current value of the progressBar + addStep(3) will add 3 steps """ - plusOne=self.prgBar.value()+1 + plusOne=self.prgBar.value()+step self.prgBar.setValue(plusOne) def reset(self): """!@brief Simply remove progressBar and reset cursor diff --git a/ui/dzetsaka_dock.py b/ui/dzetsaka_dock.py index 1be1553..1801e50 100644 --- a/ui/dzetsaka_dock.py +++ b/ui/dzetsaka_dock.py @@ -2,7 +2,7 @@ # Form implementation generated from reading ui file 'dzetsaka_dock.ui' # -# Created: Sat May 28 22:59:25 2016 +# Created: Sun Jun 12 09:26:04 2016 # by: PyQt4 UI code generator 4.10.4 # # WARNING! All changes made in this file will be lost!