From 34beacef7a9fd7ab56b98f33acefa564fffd969b Mon Sep 17 00:00:00 2001
From: Commandre Benjamin <benjamin.commandre@irstea.fr>
Date: Tue, 18 Sep 2018 18:01:53 +0200
Subject: [PATCH] Rasterisation avec OTB
---
Popup.py | 10 ++--
Processing.py | 3 +-
RasterSat_by_date.py | 113 +------------------------------------------
Vector.py | 65 ++++++++-----------------
4 files changed, 27 insertions(+), 164 deletions(-)
diff --git a/Popup.py b/Popup.py
index ea9e1a0..be53301 100644
--- a/Popup.py
+++ b/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.
diff --git a/Processing.py b/Processing.py
index edcd874..0a1b325 100644
--- a/Processing.py
+++ b/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)
diff --git a/RasterSat_by_date.py b/RasterSat_by_date.py
index e40f0fc..bb32a03 100644
--- a/RasterSat_by_date.py
+++ b/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()
-
-
- 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
+ otb_MathBand.ExecuteAndWriteOutput()
\ No newline at end of file
diff --git a/Vector.py b/Vector.py
index 2569aaf..e9fd31c 100644
--- a/Vector.py
+++ b/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
--
GitLab