v1.1

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()
 

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

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)

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()

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

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!