En la clasificación por pixel, cada pixel en la imágen es un elemento aislado y esperamos clasificarlos a partir de la información que nos da cada uno independientemente. Esta técnica es más usada cuando se tiene suficiente información temporal, es decir, se por pixel pero en un periodo de tiempo, con lo cual se pueden obtener estadísticas como máximos, valores medios de las bandas y otros índices que se obtienen a partir de las bandas (principalmente NDVI y NWDI). Sin embargo, si no se tiene buena cobertura temporal entonces este enfoque es limitado.
En este enfoque primero se hace una segmentación de la imágen y se busca clasificar los segmentos (también referidos como objetos). Esta es la técnica que vamos a usar en la tarea. Para hacer la segmentación voy a usar SLIC (http://ivrg.epfl.ch/research/superpixels) con la implementación de scikit-image (http://scikit-image.org/, http://scikit-image.org/docs/0.10.x/api/skimage.segmentation.html#slic).
Para nuestro ejemplo, vamos a usar la implementación de SLIC en scikit-image. Primero veremos un ejemplo de cómo usar SLIC con una imagen Landsat (originalmente las imágenes de Landsat vienen separadas en archivos, donde cada uno pertence a una de las bandas del sensor, en este ejemplo usamos una imagen que fue preparada previamente y que ya tiene todas las bandas en el mismo raster).
In [ ]:
%matplotlib inline
from skimage.segmentation import slic
from osgeo import gdal
import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from skimage import io
from glob import glob as g
In [ ]:
landsat_files = "../datos/ejercicio_final/rasters/limpios/ls_edo_mex_%i.tif"
raster = gdal.Open(landsat_files % 0)
num_bands = raster.RasterCount
num_rows = raster.RasterYSize
num_cols = raster.RasterXSize
#Necesitamos extraer los datos y ponerlos en un arreglo de dimensiones (ancho, largo, 3) para la función slic
im_data = np.zeros((num_rows, num_cols, num_bands), np.float)
for i in range(num_bands):
im_data[:,:,i] = raster.GetRasterBand(i+1).ReadAsArray()
In [ ]:
slic_data = im_data[2000:2250, 2000:2250]
segments = slic(slic_data, convert2lab = False, multichannel = True, compactness = 15, sigma = 0.4, n_segments = 200)
print segments.max()
In [ ]:
plt.figure(figsize = (15,15))
io.imshow(segments, cmap = cm.Set3)
In [ ]:
from skimage import io
from skimage.segmentation import mark_boundaries
Normalizamos para visualizar (algo así hace QGIS en la opción de visualización que usa valor medio y desviación estándar)
In [ ]:
bounds_color = (1, 1, 1)
means = slic_data[:,:,(3, 1, 0)].reshape((-1,3)).mean(0)
stdv = slic_data[:,:,(3, 1, 0)].reshape((-1,3)).std(0)
rgb_data = (slic_data[:,:,(3, 1, 0)] - means) / (2*stdv)
rgb_data[rgb_data > 1] = 1
rgb_data[rgb_data < -1] = -1
rgb_data -= rgb_data.reshape((-1,3)).min(0)
rgb_data /= rgb_data.reshape((-1,3)).max(0)
plt.figure(figsize = (15, 15))
#mark_boundaries marca los bordes de los segmentos
io.imshow(mark_boundaries(rgb_data, segments, bounds_color, None))
In [ ]:
from sklearn.decomposition import PCA, FastICA
In [ ]:
pca = PCA()
pca.fit(im_data.reshape(-1, num_bands))
In [ ]:
pca_data = pca.transform(slic_data.reshape(-1, num_bands))
In [ ]:
print pca_data.max()
print pca_data.min()
In [ ]:
pca_im = pca_data[:, :3].reshape(250, 250, -1)
segments = slic(pca_im, convert2lab = False, multichannel = True, compactness = 15, sigma = 0.4, n_segments = 200)
In [ ]:
segments.max()
In [ ]:
plt.figure(figsize = (15, 15))
io.imshow(segments, cmap = cm.Set3)
In [ ]:
plt.figure(figsize = (15, 15))
io.imshow(mark_boundaries(rgb_data, segments, bounds_color, None))
El ejemplo anterior, es sólo para ver cómo funciona la segmentación, a continuación procesamos todas las imágenes, las cuales fueron preparadas como se describe en el Apéndice.
In [ ]:
#Ahora generamos las segmentaciones
num_total_segs = 0
recortes = glob(con_etiquetas_dir + "/*.tif")
print con_etiquetas_dir
print len(recortes)
for recorte in recortes:
data_src = gdal.Open(recorte)
xsize = data_src.RasterXSize
ysize = data_src.RasterYSize
num_bands = data_src.RasterCount
data = np.zeros((ysize, xsize, num_bands), np.float)
for i in range(num_bands):
data[:,:,i] = data_src.GetRasterBand(i+1).ReadAsArray()
segments = slic(data, convert2lab = False, multichannel = True, compactness = 15, sigma = 0.4, n_segments = 600)
segs_src = data_src.GetDriver().Create(recorte.replace(".tif", ".segs"), xsize, ysize, 1, gdal.GDT_Int16)
segs_src.SetProjection(data_src.GetProjection())
segs_src.SetGeoTransform(data_src.GetGeoTransform())
segs_src.GetRasterBand(1).WriteArray(segments)
segs_src = None
data_src = None
num_total_segs += len(np.unique(segments))
print "Promedio segmentos por imagen: %f" % (num_total_segs / len(recortes))
Para validar
Vamos a trabajar en el Estado de México.
In [ ]:
import glob as g
from osgeo import gdal
def unir_bandas_ls(src_dir, dst_name):
band_files = g.glob(src_dir + "*_B*.TIF")[:-1]
# no vamos a usar la banda 8 (pancromático) porque tiene una resolución más alta
print '\n'.join(band_files)
#inciar nuevo raster con los datos de la primer banda de la imagen landsat que queremos unir
raster_band = gdal.Open(band_files[0])
driver = raster_band.GetDriver()
dst_raster = driver.Create(dst_name, raster_band.RasterXSize, raster_band.RasterYSize, 8, gdal.GDT_Byte)
dst_raster.SetProjection(raster_band.GetProjection())
dst_raster.SetGeoTransform(raster_band.GetGeoTransform())
for i, f in enumerate(band_files):
raster_band = gdal.Open(f)
band_data = raster_band.GetRasterBand(1).ReadAsArray()
dst_raster.GetRasterBand(i+1).WriteArray(band_data)
dst_raster = None
In [ ]:
def gdal_unir_bandas_ls(src_dir, dst_name):
band_files = g.glob(src_dir + "*_B*.TIF")[:-1]
# no vamos a usar la banda 8 (pancromático) porque tiene una resolución más alta
!gdalbuildvrt -separate {dst_name.replace('tif', 'vrt')} {' '.join(band_files)}
!gdal_translate {dst_name.replace('tif', 'vrt')} $dst_name
!rm {dst_name.replace('tif', 'vrt')}
In [ ]:
src_dirs = ["../datos/landsat/edo_mexico/934565/2000/2000-01-24/l1t/",
"../datos/landsat/edo_mexico/934566/2000/2000-01-24/l1t/",
"../datos/landsat/edo_mexico/967057/2000/2000-01-01/l1t/",
"../datos/landsat/edo_mexico/967058/2000/2000-02-02/l1t/"]
dst_name = "../datos/ejercicio_final/rasters/ls_edo_mex_%i.tif"
for i, src_dir in enumerate(src_dirs):
print "Create: " + dst_name % i
gdal_unir_bandas_ls(src_dir, dst_name % i)
In [ ]:
# Recortar imagenes para remover las orillas que no sirven. Para esto usamos el mosaico del mundo http://landsat.usgs.gov/documents/wrs2_descending.zip
cut_shp = "../datos/ejercicio_final/tiles_shapes/ls_edo_mex_%i.shp" #estos son shapes que corresponden a cada una de las zonas landsat, fueron extraídos manualmente
src_tif = "../datos/ejercicio_final/rasters/ls_edo_mex_%i.tif"
dst_tif = "../datos/ejercicio_final/rasters/limpios/ls_edo_mex_%i.tif"
for i in range(4):
print "Warp %i" % i
!gdalwarp -dstnodata 0 -q -cutline {cut_shp % i} -crop_to_cutline -of GTiff {src_tif % i} {dst_tif % i}
In [ ]:
src_file = "../datos/ejercicio_final/ls_edo_mex_%i.tif"
!gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o ../datos/ejercicio_final/prueba_merge.tif {src_file % 0} {src_file % 1} {src_file % 2} {src_file % 3}
Generamos un shape con números en vez de caracteres para las clases, para poder escribir los valores en un raster
In [ ]:
import fiona
orig_shape_file = "../datos/ejercicio_final/labels/edo_mex_labels.shp"
int_shape_file = "../datos/ejercicio_final/labels/edo_mex_int_labels.shp"
src_tif = "../datos/ejercicio_final/rasters/limpios/ls_edo_mex_%i.tif"
In [ ]:
with fiona.open(orig_shape_file) as src:
src_driver = src.driver
src_crs = src.crs
src_schema = src.schema
recs = list(src)
print "El esquema actual: " + str(src_schema['properties'])
src_schema['properties']['level_1'] = "int:8"
print "El nuevo esquema: " + str(src_schema['properties'])
#Queremos que el nivel 1 sea numérico en vez de una cadena de caracteres
with fiona.open(int_shape_file,
'w',
driver=src_driver,
crs=src_crs,
schema=src_schema) as c:
clases_l1 = {u'TIERRAS FORESTALES': 1,
u'INDEFINIDO': 8,
u'PRADERAS': 3,
u'TIERRAS DE USO AGRICOLA': 4,
u'AGUA': 5,
u'HUMEDAL': 2,
u'ASENTAMIENTOS': 6,
u'OTROS': 7
}
for rec in recs:
clase = rec["properties"]["level_1"]
if clase not in clases_l1:
clases_l1[clase] = len(clases_l1.keys())
rec["properties"]["level_1"] = clases_l1[clase]
c.write(rec)
print clases_l1
c.closed
Ahora vamos a escribir las etiquetas y los ids de segmentos de la clasificación a archivos raster. Para esto:
In [ ]:
format = "GTiff"
driver = gdal.GetDriverByName( format )
oids_tif = "../datos/ejercicio_final/rasters/limpios/ls_edo_mex_%i.oids"
labels_tif = "../datos/ejercicio_final/rasters/limpios/ls_edo_mex_%i.gt"
labels_src = int_shape_file
for i in range(4):
print src_tif % i
ref_rast = gdal.Open(src_tif % i) #raster de referencia para crear el raster destino
blank_rast = driver.Create(oids_tif % i, ref_rast.RasterXSize, ref_rast.RasterYSize, 1, gdal.GDT_Int16)
blank_rast.SetProjection(ref_rast.GetProjection())
blank_rast.SetGeoTransform(ref_rast.GetGeoTransform())
blank_rast = None
for i in range(4):
!gdal_rasterize -a OBJECTID -l edo_mex_int_labels $labels_src {oids_tif % i}
Ahora creamos rasters con las etiquetas de clase.
In [ ]:
for i in range(4):
ref_rast = gdal.Open(src_tif % i)
blank_rast = driver.Create(labels_tif % i, ref_rast.RasterXSize, ref_rast.RasterYSize, 1, gdal.GDT_Byte)
blank_rast.SetProjection(ref_rast.GetProjection())
blank_rast.SetGeoTransform(ref_rast.GetGeoTransform())
blank_rast = None
for i in range(4):
!gdal_rasterize -a level_1 -l edo_mex_int_labels $labels_src {labels_tif % i}
Ahora subdividimos las imágenes en imágenes de 500x500 pixeles
In [ ]:
target_dir = "../datos/ejercicio_final/recortes/500x500/"
src_file = "../datos/ejercicio_final/rasters/limpios/ls_edo_mex_%i.%s"
tile_size = 500
for i in range(4):
for ext in ('tif', 'gt', 'oids'):
if not os.path.exists(target_dir + ext):
os.mkdir(target_dir + ext)
!gdal_retile.py -targetDir {target_dir + ext} -ps $tile_size $tile_size {src_file % (i, ext)}
out_files = g.glob(target_dir + ext + "/*.tif")
for filename in out_files:
newname = filename.replace(ext + "/", "").replace(".tif", "." + ext)
os.rename(filename, newname)
for ext in ('tif', 'gt', 'oids'):
os.rmdir(target_dir + ext)
Depuramos el conjunto de imágenes para confirmar que no contengan pixeles sin datos (ceros) y que todas tengan al menos 100 pixeles anotados.
In [ ]:
import os
recortes = g.glob("../datos/ejercicio_final/recortes/500x500/*.tif")
num_nodata = 0
sin_etiqueta = 0
for recorte in recortes:
r = gdal.Open(recorte)
l = gdal.Open(recorte.replace(".tif",".gt"))
ldata = l.GetRasterBand(1).ReadAsArray()
rdir = os.path.dirname(recorte)
num_bands = r.RasterCount
nrows = r.RasterYSize
ncols = r.RasterXSize
nodata_percent = .8
has_data = True
labels = False
for i in range(num_bands):
bdata = r.GetRasterBand(i+1).ReadAsArray()
if (bdata == 0).sum() >= 1:
has_data = False
break
if has_data:
if ldata.sum() > 100:
labels = True
else:
sin_etiqueta += 1
else:
num_nodata += 1
# sin_etiquetas_dir = os.path.join(rdir, "sin_etiquetas")
# if not os.path.exists(sin_etiquetas_dir):
# !mkdir $sin_etiquetas_dir
# !cp $recorte {os.path.join(sin_etiquetas_dir, os.path.basename(recorte))}
if has_data and labels:
con_etiquetas_dir = os.path.join(rdir, "con_etiquetas_100")
if not os.path.exists(con_etiquetas_dir):
os.mkdir(con_etiquetas_dir)
!cp $recorte {os.path.join(con_etiquetas_dir, os.path.basename(recorte))}
labels_file = recorte.replace(".tif",".gt")
!cp $labels_file {os.path.join(con_etiquetas_dir, os.path.basename(labels_file))}
oids_file = recorte.replace(".tif",".oids")
!cp $oids_file {os.path.join(con_etiquetas_dir, os.path.basename(oids_file))}
print "%i archivos sin datos" % num_nodata
print "Sin etiquetas %i" % sin_etiqueta
Checamos el set de datos.
Validación
In [ ]:
def kappa(conf_mat):
pr_e = (conf_mat.sum(1)*conf_mat.sum(0)).sum()/float(conf_mat.sum())**2
pr_a = conf_mat.diagonal().sum()/float(conf_mat.sum())
return (pr_a - pr_e) / (1 - pr_e)
def accuracy(conf_mat):
return conf_mat.diagonal().sum() / conf_mat.astype(float).sum()
def precision(conf_mat):
return conf_mat.diagonal() / conf_mat.astype(float).sum(1)
def cobertura(conf_mat):
return conf_mat.diagonal() / conf_mat.astype(float).sum(0)
def f1(conf_mat):
return (precision(conf_mat)*cobertura(conf_mat)) / (precision(conf_mat) + cobertura(conf_mat))