Clasificación

Por pixel

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.

Por objeto o segmento

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).

Segmentación de imágenes

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))

Validación

Para validar

Referencias

Apéndice: Lectura de imágenes landsat y generación de etiquetas en raster a partir del shapefile

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}

Pegar rasters en mosaico


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}

Escribir etiquetas y ids de objetos interpretados manualmente en rasters

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:

  1. Generamos primero un raster sin datos que tenga las dimensiones y la proyección que corresponde a cada zona.
  2. Rasterizamos el shape de etiquetas escribiendo al raster que recien creamos. De este modo escribimos sólo los datos que corresponden a la zona que abarca un tile. En este caso hay cuatro tiles que cubren alguna parte del Estado de México.

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))