Processing.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# This file is part of PHYMOBAT 2.0.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
# 
# PHYMOBAT 2.0 is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# PHYMOBAT 2.0 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with PHYMOBAT 2.0.  If not, see <http://www.gnu.org/licenses/>.

import os, sys
import numpy as np
import configparser
import glob 
from osgeo import gdal, ogr, osr

import otbApplication as otb

# Class group image
from app.Archive import Archive

import app.Constantes as Constantes
import app.Outils as Outils

from collections import defaultdict
from multiprocessing import Process
from multiprocessing.managers import BaseManager, DictProxy

class Processing(object):
	
	def __init__(self):
		# List of output raster path
		self.raster_path = defaultdict(dict)
				
	def i_download(self):
		"""
			Méthode pour télécharger les images sur le site Theia Land.
		"""

		self.images_folder = "{0}/Images".format(self.path_folder_dpt)

		if not os.path.exists(self.images_folder) :
			os.makedirs(self.images_folder)

		self.check_download = Archive(self.captor_project, self.classif_year, self.path_area, self.images_folder, self.object_area)
		self.nb_avalaible_images = self.check_download.listing()
		self.check_download.download_auto(self.user, self.password, self.proxy)

		self.liste_dossier = dict()

		for dossier in sorted(glob.glob("{0}/*".format(self.images_folder))) :
			if os.path.isdir(dossier) :
				self.liste_dossier[os.path.basename(dossier)] = sorted([ x for x in glob.glob("{0}/*".format(dossier)) if x.endswith(".tif")] )
	

	def polygonisation(self, dossier, image_seuil):
		"""
			Méthode effactuant la polygonisation des images seuillées
		"""	
		
		dossier_resultats = "{0}/Resultats/{1}".format(self.path_folder_dpt, dossier)
		sortie = "{0}/{1}.shp".format(dossier_resultats, os.path.basename(image_seuil)[:-4])

		if not os.path.exists(dossier_resultats) :
			os.makedirs(dossier_resultats)

		format			= "ESRI Shapefile"
		dst_layername 	= None
		dst_fieldname 	= "Type"
		dst_field 		= 0
		src_band_n 		= 1

		src_ds = gdal.Open(image_seuil)
		srcband = src_ds.GetRasterBand(src_band_n)

		maskband = srcband.GetMaskBand()

		drv = ogr.GetDriverByName(format)
		dst_ds = drv.CreateDataSource(sortie)

		srs = osr.SpatialReference()
		srs.ImportFromWkt(src_ds.GetProjectionRef())

		dst_layer = dst_ds.CreateLayer('out', geom_type=ogr.wkbPolygon, srs=srs)

		fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
		dst_layer.CreateField(fd)
		
		gdal.Polygonize( srcband, maskband, dst_layer, dst_field)

	def i_images_processing(self): 
		"""
			Méthode appelant le calcul du MSAVI2 via OTB
			et effectuant le seuillage via MathBandX (OTB) 
		"""

		dossier_etude = "{0}/Zone_etude".format(self.path_folder_dpt)
		dossier_MSAVI = "{0}/MSAVI2".format(self.path_folder_dpt)
		dossier_seuil = "{0}/Seuillage".format(self.path_folder_dpt)

		if not os.path.exists(dossier_etude) :
			os.makedirs(dossier_etude)

		if not os.path.exists(dossier_MSAVI) :
			os.makedirs(dossier_MSAVI)

		if not os.path.exists(dossier_seuil) :
			os.makedirs(dossier_seuil)
		
		otb_MathBand = otb.Registry.CreateApplication("BandMathX")
		otb_MathBand.SetParameterString("exp", "(im1b1 > -10000 ? (im1b1 > 0.025 ? (im1b1 < 0.15 ? 2 : 3 ) : 1) : im1b1)" )

		otb_MSAVI = otb.Registry.CreateApplication("RadiometricIndices")
		otb_MSAVI.SetParameterInt("channels.blue", 1)
		otb_MSAVI.SetParameterInt("channels.green", 2)
		otb_MSAVI.SetParameterInt("channels.red", 3)
		otb_MSAVI.SetParameterInt("channels.nir", 4)
		otb_MSAVI.SetParameterStringList("list", ["Vegetation:MSAVI2"])

		for dossier in self.liste_dossier :

			nom_dst = "{0}/{1}".format(dossier_etude, dossier)
			nom_msavi = "{0}/{1}".format(dossier_MSAVI, dossier)
			nom_seuil = "{0}/{1}".format(dossier_seuil, dossier)

			if not os.path.exists(nom_dst):
				os.makedirs(nom_dst)

			if not os.path.exists(nom_msavi):
				os.makedirs(nom_msavi)

			if not os.path.exists(nom_seuil):
				os.makedirs(nom_seuil)

			for img in self.liste_dossier[dossier] :
				chemin_dst = "{0}/{1}".format(nom_dst, os.path.basename(img))
				Outils.clip(img, self.object_area, form="GTiff", dst=chemin_dst)
				
				chemin_msavi = "{0}/{1}".format(nom_msavi, os.path.basename(img))
				chemin_seuil = "{0}/{1}".format(nom_seuil, os.path.basename(img))

				otb_MSAVI.SetParameterString("in", chemin_dst)
				otb_MSAVI.SetParameterString("out", chemin_msavi)
				otb_MSAVI.ExecuteAndWriteOutput()
				# otb_MSAVI.Execute()

				otb_MathBand.SetParameterStringList("il", [chemin_msavi])
				# otb_MathBand.AddImageToParameterInputImageList("il",otb_MSAVI.GetParameterOutputImage("out"));
				otb_MathBand.SetParameterString("out", chemin_seuil)

				otb_MathBand.ExecuteAndWriteOutput()

				self.polygonisation(dossier, chemin_seuil)