#!/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
import subprocess
import numpy as np
try :
import ogr, gdal
except :
from osgeo import ogr, gdal
from rasterstats import *
from collections import *
from RasterSat_by_date import RasterSat_by_date
[docs]class Vector():
"""
Vector class to extract a area, vector data and zonal statistic (``rasterstats 0.3.2 package``)
:param vector_used: Input/Output shapefile to clip (path)
:type vector_used: str
:param vector_cut: Area shapefile (path)
:type vector_cut: str
:param data_source: Input shapefile information
:type data_source: ogr pointer
:param stats_dict: ``Rasterstats`` 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
:opt: **Remove** (int) - For the remove_shp variable
"""
def __init__(self, used, cut, **opt):
"""Create a new 'Vector' instance
"""
self.vector_cut = cut
self.vector_used = used
self.remove_shp = opt['Remove'] if opt.get('Remove') else 0
self.clip_vector()
self.data_source = ''
if self.data_source == '':
self.vector_data()
# List of field name
self.field_names = [self.data_source.GetLayer().GetLayerDefn().GetFieldDefn(l).GetName() \
for l in range(self.data_source.GetLayer().GetLayerDefn().GetFieldCount())]
self.stats_dict = defaultdict(list)
self.raster_ds = None
[docs] def clip_vector(self):
"""
Function to clip a vector with a vector
"""
outclip = os.path.split(self.vector_used)[0] + '/Clip_' + os.path.split(self.vector_used)[1]
if not os.path.exists(outclip) or self.remove_shp == 1:
print 'Clip of ' + os.path.split(self.vector_used)[1]
# Command to clip a vector with a shapefile by OGR
process_tocall_clip = ['ogr2ogr', '-overwrite', '-skipfailures', outclip, self.vector_used, '-clipsrc', self.vector_cut]
subprocess.call(process_tocall_clip)
# Replace input filename by output filename
self.vector_used = outclip
[docs] def vector_data(self):
"""
Function to extract vector layer information
"""
# import ogr variable
self.data_source = ogr.GetDriverByName('ESRI Shapefile').Open(self.vector_used, 0)
if self.data_source is None:
print('Could not open file')
sys.exit(1)
print('Shapefile opening : ' + self.data_source.GetLayer().GetName())
[docs] def close_data(self):
"""
Function to remove allocate memory
"""
# Close data sources
self.data_source.Destroy()
print('Shapefile closing : ' + self.data_source.GetLayer().GetName())
[docs] def zonal_stats(self, (inraster, band), **kwargs):
"""
Function to compute the average in every polygons for a raster
because of package ``rasterstats`` in */usr/local/lib/python2.7/dist-packages/rasterstats-0.3.2-py2.7.egg/rasterstats/*
:param (inraster,band): inraster -> Input image path, and band -> band number
:type (inraster,band): tuple
:kwargs: **rank** (int) - Zonal stats ranking launch
**nb_img** (int) - Number images launched with zonal stats
"""
ranking = kwargs['rank'] if kwargs.get('rank') else 0
nb_img = kwargs['nb_img'] if kwargs.get('nb_img') else 1
print('Compute ' + os.path.split(str(self.vector_used))[1] + ' stats on ' + os.path.split(str(inraster))[1])
stats = raster_stats(str(self.vector_used), str(inraster), stats =['mean'], band_num=band)
for i in range(len(stats)):
temp = defaultdict(lambda : [0]*nb_img)
for j in range(nb_img):
try :
temp[0][j] = self.stats_dict[i][j]
except IndexError:
pass
temp[0][ranking] = stats[i].values()[1]
self.stats_dict[i] = temp[0]
print('End of stats on ' + os.path.split(str(inraster))[1])
[docs] def zonal_stats_pp(self, inraster):
"""
A zonal statistics ++ to dertermine pxl percent in every polygon
:param inraster: Input image path
:type inraster: str
:returns: dict -- **p_stats** : dictionnary with pxl percent in every polygon. Mainly 'Maj_count' (Majority Value) and 'Maj_count_perc' (Majority Percent)
"""
p_stats = raster_stats(str(self.vector_used), str(inraster), stats=['count'], copy_properties=True, categorical=True)
for i in range(len(p_stats)):
percent = 0.0
for p in p_stats[i].keys():
if type(p) == np.float32 and p_stats[i][p]/float(p_stats[i]['count'])*100 > percent:
p_stats[i]['Maj_count'] = p
p_stats[i]['Maj_count_perc'] = p_stats[i][p]/float(p_stats[i]['count'])*100
percent = p_stats[i]['Maj_count_perc']
if not 'Maj_count' in p_stats[i].keys() or not 'Maj_count_perc' in p_stats[i].keys():
p_stats[i]['Maj_count']=0
p_stats[i]['Maj_count_perc']=0
return p_stats
[docs] 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.
: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
: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('', '', [0]) # Call the raster class
example_raster.choice_nb_b = kwargs['choice_nb_b'] if kwargs.get('choice_nb_b') else 0
raster_info = example_raster.raster_data(raster_head)# Extract data info
# Define the validation's vector
valid_raster = self.vector_used[:-3]+'TIF' # Name of the output raster
if os.path.exists(str(valid_raster)):
os.remove(valid_raster)
# Create the empty raster with the same properties
# Condition for the rasters with several bands
if raster_info[1].RasterCount > 1:
data_raster = raster_info[0][0]
else:
data_raster = raster_info[0]
info_out = example_raster.create_raster(valid_raster, data_raster, raster_info[1])
self.raster_ds = example_raster.out_ds
# Virtual rasterize the vector
pt_rast = gdal.RasterizeLayer(self.raster_ds, [1], self.data_source.GetLayer(), \
options=["ATTRIBUTE=" + str(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)
return valid_raster