Rasterisation avec OTB

Popup.py

@@ -83,10 +83,12 @@ class proxy_window(QtWidgets.QWidget):
 		"""
 			Function to use input proxy id
 		"""
-		self.login_proxy = "%s" % self.w_proxy.lineEdit_login_proxy.text()
-		self.password_proxy = "%s" % self.w_proxy.lineEdit_password_proxy.text()
-		self.proxy = "%s" % self.w_proxy.lineEdit_proxy.text()
-	
+		self.login_proxy = "{0}".format(self.w_proxy.lineEdit_login_proxy.text())
+		self.password_proxy = "{0}".format(self.w_proxy.lineEdit_password_proxy.text())
+		self.proxy = "{0}".format(self.w_proxy.lineEdit_proxy.text())
+
+		self.close_window()
+
 	def close_window(self):
 		"""
 			Function to close the popup.

Processing.py

@@ -709,9 +709,8 @@ class Processing(object):
 		rpg_tif.clip_vector(os.path.dirname(self.sample_name[0]))
 		rpg_tif.vector_data()
 
-		kwargs['choice_nb_b'] = 1
 		self.logger.debug(self.path_ortho)
-		out_carto.stats_rpg_tif = out_carto.zonal_stats_pp(rpg_tif.layer_rasterization(self.path_ortho, 'CODE_GROUP', **kwargs))
+		out_carto.stats_rpg_tif = out_carto.zonal_stats_pp(rpg_tif.layer_rasterization(self.path_ortho, 'CODE_GROUP'))
 
 		# Final cartography
 		out_carto.create_cartography(self.out_fieldname_carto, self.out_fieldtype_carto)

RasterSat_by_date.py

@@ -306,115 +306,4 @@ class RasterSat_by_date(object):
 		.format(canal_PIR, canal_R))
 		otb_MathBand.SetParameterString("out", img_ndvi)
 
-		otb_MathBand.ExecuteAndWriteOutput()
-
+		otb_MathBand.ExecuteAndWriteOutput()
\ No newline at end of file
-
-	def create_empty_raster(self, out_raster_path, in_ds):
-		"""
-			Create an empty raster with the input raster property
-			
-			@param out_raster_path: Output image path
-			@type out_raster_path: str
-
-			@param in_ds: Raster information
-			@type in_ds: gdal pointer
-			
-			@returns: gdal pointer -- variable **out_ds**, Raster out information.
-					  
-					  int -- variable **nbband**, Band number of the out layer. 
-					  
-					  int -- variable **e**, Index to know if the raster exists.
-					  If it doesn't exists e = 0 else e = 1 (by default).
-		"""
-		
-		e = 1 # Raster out exists by default 
-		
-		# Verify if the processing take input band or one spectral band	
-		if self.choice_nb_b == 1:
-			nbband = 1
-		else:
-			nbband = in_ds.RasterCount 
-			
-		driver = osgeo.gdal.GetDriverByName('GTiff') 
-
-		if os.path.exists(out_raster_path):
-			os.remove(out_raster_path)
-
-		if not os.path.exists(out_raster_path):
-
-			e = 0	
-
-			# Create outfile
-			self.out_ds = driver.Create(out_raster_path, in_ds.RasterXSize, in_ds.RasterYSize, nbband, osgeo.gdal.GDT_Float32)
-
-			if self.out_ds is None:
-				self.logger.error('Could not create {0}'.format(os.path.basename(out_raster_path)))
-				sys.exit(1)
-				
-			# Get extent coordinates and raster resolution
-			transform = in_ds.GetGeoTransform()
-
-			minX = transform[0]
-			maxY = transform[3]
-			pixelWidth = transform[1]
-			pixelHeight = transform[5]
-			
-			geotransform = [minX, pixelWidth, 0, maxY, 0, pixelHeight]
-			
-			# Record projection
-			def_projection = in_ds.GetProjection() 
-
-			# Set the geo-traking and outfile projection
-			self.out_ds.SetGeoTransform(geotransform)
-			self.out_ds.SetProjection(def_projection)
-
-		else:
-			self.out_ds = osgeo.gdal.Open(out_raster_path, osgeo.gdal.GA_ReadOnly)
-			
-		return nbband, e
-	
-	def complete_raster(self, nbband, e, data): 
-		"""
-			This function complete the function above : func:`create_empty_raster()`. 
-			It fills the raster table and close the layer.
-		
-			@param out_ds: Raster out information
-			@type out_ds: gdal pointer
-
-			@param nbband: Band number of the out layer
-			@type nbband: int
-
-			@param e: Index to know if the raster existed. If it didn't exist e = 0.
-			@type e: int
-
-			@param data: Pixel value matrix. Matrix size equal to that of a raster.
-			@type data: numpy.array
-		"""
-		
-		# The e index to verify if the layer existed already because of the 
-		# function :func:`create_empty_raster()`
-		if e == 0 :   
-			
-			p = 0 # Increment for the number band
-			while p < nbband:
-				#Incrementation
-				p = p + 1
-			
-				self.logger.info("Copy on the band : {0}".format(p))
-	  
-				# Loading spectral band of outfile
-				out_band = self.out_ds.GetRasterBand(p) 
-				# write the data
-				if data.ndim == 2:
-					out_band.WriteArray(data, 0, 0)
-				else:
-					out_band.WriteArray(data[p-1], 0, 0)
-				
-				# Closing and statistics on output raster
-				out_band.FlushCache()
-				# out_band.SetNoDataValue(-10000)
-				out_band.GetStatistics(-1, 1) 
-				out_band = None	
-			
-		# Close open data
-		self.out_ds = None
\ No newline at end of file

Vector.py

@@ -21,14 +21,12 @@ import os, sys
 import subprocess
 import numpy as np
 import osgeo
-import time
 
 # try :
 # 	import ogr, gdal
 # except :
 # 	from osgeo import ogr, gdal
 	
-import rasterstats
 import otbApplication as otb
 from collections import *
 import Outils
@@ -36,7 +34,7 @@ from RasterSat_by_date import RasterSat_by_date
 
 class Vector():
 	"""	
-		Vector class to extract a area, vector data and zonal statistic (``rasterstats 0.3.2 package``)
+		Vector class to extract a area, vector data and zonal statistic
 		
 		@param vector_used: Input/Output shapefile to clip (path)
 		@type vector_used: str
@@ -53,12 +51,9 @@ class Vector():
 		@param data_source: Input shapefile information
 		@type data_source: ogr pointer
 
-		@param stats_dict: ``Rasterstats`` results
+		@param stats_dict: Stats results
 		@type stats_dict: dict
 
-		@param raster_ds: Raster information for a probably rasterization
-		@type raster_ds: gdal pointer
-
 		@param remove_shp: Remove shapefile or not. 0 : don't remove, 1 : remove
 		@type remove_shp: int
 
@@ -81,7 +76,6 @@ class Vector():
 		self.remove_shp = opt['Remove'] if opt.get('Remove') else 0
 	
 		self.stats_dict = defaultdict(list)
-		self.raster_ds = None
 	
 	def clip_vector(self, output_folder):
 		"""
@@ -202,53 +196,32 @@ class Vector():
 				
 	def layer_rasterization(self, raster_head, attribute_r, **kwargs):
 		"""
-			Function to rasterize a vector. 
-			Define resolution, projection of the output raster with a raster head.
-			And complete the gdal pointer empty properties with the layer's information
-			 of the vector and a defined field.
-			If a raster has several band, in option you can choice if you want one band or more. 
+			Function to rasterize a vector using OTB.
 			
-			:param raster_head: Raster path that will look like the final raster of the rasterization
-			:type raster_head: str
+			@param raster_head: Raster path that will look like the final raster of the rasterization
+			@type raster_head: str
 
-			:param attribute_r: Field name of the shapefile that contains class names
-			:type attribute_r: str
+			@param attribute_r: Field name of the shapefile that contains class names
+			@type attribute_r: str
 
-			:kwargs: **choice_nb_b** (int) - Output image number of band. If you choice 1, take first band. If you choice 2, take two first band etc... 
-			:returns: str -- **valid_raster** : output raster path from the rasterization
+			@kwargs: **choice_nb_b** (int) - Output image number of band. If you choice 1, take first band. If you choice 2, take two first band etc... 
+			@returns: str -- **valid_raster** : output raster path from the rasterization
 		"""
-		
-		# Export a example of a raster out information
-		# for the validation shapefile
-		example_raster = RasterSat_by_date('', '') # Call the raster class
-		example_raster.choice_nb_b = kwargs['choice_nb_b'] if kwargs.get('choice_nb_b') else 0
 
-		in_ds = osgeo.gdal.Open(raster_head, osgeo.gdal.GA_ReadOnly)
-		
-		# Define the validation's vector
 		valid_raster = "{0}.TIF".format(self.vector_used[:-4])# Name of the output raster
+
+		self.logger.debug("valid_raster : {0}".format(valid_raster))
 		
 		if os.path.exists(valid_raster):
 			os.remove(valid_raster)
-			
-		info_out = example_raster.create_empty_raster(valid_raster, in_ds)
-
-		in_ds = None
-
-		self.raster_ds = example_raster.out_ds
-		
-		# Virtual rasterize the vector 
-		pt_rast = osgeo.gdal.RasterizeLayer(self.raster_ds, [1], self.data_source.GetLayer(), \
-									  options=["ATTRIBUTE = {0}".format(attribute_r)])
-		
-		if pt_rast != 0:
-			raise Exception("error rasterizing layer: %s" % pt_rast)
 		
-		new_data = self.raster_ds.ReadAsArray()
-		
-		self.raster_ds = None
-		# Complete the raster creation
-		example_raster.complete_raster(*info_out, new_data)
-		example_raster.logger.close()
+		layerRasterization = otb.Registry.CreateApplication("Rasterization")
+		layerRasterization.SetParameterString("in", self.vector_used)
+		layerRasterization.SetParameterString("out", valid_raster)
+		layerRasterization.SetParameterString("im", raster_head)
+		layerRasterization.SetParameterString("mode", "attribute")
+		layerRasterization.SetParameterString("mode.attribute.field", attribute_r)
+
+		layerRasterization.ExecuteAndWriteOutput()
 		
 		return valid_raster
\ No newline at end of file