Code source de Toolbox
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# This file is part of PHYMOBAT 1.2.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
#
# PHYMOBAT 1.2 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 1.2 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 1.2. If not, see <http://www.gnu.org/licenses/>.
import os, sys, subprocess, glob
import numpy as np
from datetime import date
import json
[docs]class Toolbox():
"""
Class used to grouped small tools to cut raster or compute statistics on a raster.
:param imag: Input image (path)
:type imag: str
:param vect: Extent shapefile (path)
:type vect: str
"""
def __init__(self):
"""
Create a new 'Toolbox' instance
"""
self.imag = ''
self.vect = ''
[docs] def clip_raster(self, **kwargs):
"""
Function to clip a raster with a vector. The raster created will be in the same folder than the input raster.
With a prefix *Clip_*.
:kwargs: **rm_rast** (int) - 0 (by default) or 1. Variable to remove the output raster. 0 to keep and 1 to remove.
:returns: str -- variable **outclip**, output raster clip (path).
"""
rm_rast = kwargs['rm_rast'] if kwargs.get('rm_rast') else 0
outclip = os.path.split(str(self.imag))[0] + '/Clip_' + os.path.split(str(self.imag))[1]
if not os.path.exists(outclip) or rm_rast == 1:
print 'Raster clip of ' + os.path.split(str(self.imag))[1]
# Command to clip a raster with a shapefile by Gdal
process_tocall_clip = ['gdalwarp', '-overwrite', '-dstnodata', '-10000', '-q', '-cutline', self.vect, '-crop_to_cutline', '-of', 'GTiff', self.imag, outclip]
# This is a trick to remove warning with the polygons that touch themselves
try:
r = subprocess.call(process_tocall_clip)
if r == 1:
sys.exit(1)
except SystemExit:
print("Dissolve vector cut of the validation !!!")
vect_2 = self.vect[:-4] + '_v2.shp'
preprocess_tocall = 'ogr2ogr -overwrite ' + vect_2 + ' ' + self.vect + ' -dialect sqlite -sql "SELECT ST_Union(geometry), * FROM ' + \
os.path.basename(self.vect)[:-4] +'"'
os.system(preprocess_tocall)
print 'Raster clip of ' + os.path.split(str(self.imag))[1]
process_tocall_clip = ['gdalwarp', '-overwrite', '-dstnodata', '-10000', '-q', '-cutline', vect_2, '-crop_to_cutline', '-of', 'GTiff', self.imag, outclip]
subprocess.call(process_tocall_clip)
for rem in glob.glob(vect_2[:-4] + '*'):
os.remove(rem)
return outclip
[docs] def calc_serie_stats(self, table):
"""
Function to compute stats on temporal cloud and ndvi spectral table
Ndvi stats : min max std max-min
:param table: Spectral data, cloud raster and ndvi raster
:type table: numpy.ndarray
:returns: list of numpy.ndarray -- variable **account_stats**, list of temporal NDVI stats.
numpy.ndarray -- variable **account_cloud**, pixel number clear on the area.
"""
# Compute stats on these indexes
ind = ['np.min(tab_ndvi_masked, axis=2)', 'np.max(tab_ndvi_masked, axis=2)', 'np.std(tab_ndvi_masked, axis=2)', \
'np.max(tab_ndvi_masked, axis=2)-np.min(tab_ndvi_masked, axis=2)'] # [Min, Max, Std, Max-Min]
# For the cloud map
# In the input table the cloud floor is the 5th
tab_cloud = np.dstack(table[5]) # Stack cloud table (dimension : 12*X*Y to X*Y*12)
cloud_true = (tab_cloud == 0) # if tab_cloud = 0 then True else False / Mask cloud
account_cloud = np.sum(cloud_true, axis=2) # Account to tab_cloud if != 0 => Sum of True. (Dimension X*Y)
# For the ndvi stats
# In the input table the ndvi floor is the 7th
stack_ndvi = np.dstack(table[7]) # Like cloud table, stack ndvi table
# mask_ndvi = np.ma.masked_equal(stack_ndvi, -10000, copy=True) # Mask values -10000
mask_cloud = np.ma.masked_where((stack_ndvi == -10000) | (tab_cloud != 0), tab_cloud) # Mask values -10000 and > 0(Not cloud)
# mask_cloud = (tab_cloud != 0) | (stack_ndvi == -10000)
tab_ndvi_masked = np.ma.array(stack_ndvi, mask=mask_cloud.mask) # mask_cloud.mask) # Ndvi table with clear values
# Stats on the indexes defined above
account_stats = []
for i in ind:
i_stats = eval(i) # Compute stats
i_stats.fill_value = -10000 # Substitute default fill value by -10000
account_stats.append(i_stats.filled()) # Add stats table with true fill value
# To extract index of the minimum ndvi to find the date
if ind.index(i) == 0:
i_date = eval(i.replace('np.min','np.argmin'))
# To create a matrix with the date of pxl with a ndvi min
tab_date = np.transpose(np.array(table[:3], dtype=object))
# Loop on the temporal sequency
for d in range(len(tab_date)):
mask_data = np.ma.masked_values(i_date, d)
i_date = mask_data.filled(date(int(tab_date[d][0]), int(tab_date[d][1]), int(tab_date[d][2])).toordinal()) # Date = day for year =1 day = 1 and month = 1
account_stats.append(i_date) # Add the date table
return account_stats, account_cloud
[docs] def check_proj(self):
"""
Function to check if raster's projection is RFG93.
For the moment, PHYMOBAT works with one projection only Lambert 93 EPSG:2154
"""
# Projection for PHYMOBAT
epsg_phymobat = '2154'
proj_phymobat = 'AUTHORITY["EPSG","' + epsg_phymobat + '"]'
info_gdal = 'gdalinfo -json ' + self.imag
info_json = json.loads(subprocess.check_output(info_gdal, shell=True))
# Check if projection is in Lambert 93
if not proj_phymobat in info_json['coordinateSystem']['wkt']:
output_proj = self.imag[:-4] + '_L93.tif'
reproj = 'gdalwarp -t_srs EPSG:' + epsg_phymobat + ' ' + self.imag + ' ' + output_proj
os.system(reproj)
# Remove old file and rename new file like the old file
os.remove(self.imag)
os.rename(output_proj, self.imag)