Configuración inicial:
In [1]:
%matplotlib inline
from osgeo import gdal
import gdalnumeric
import osgeo.gdalconst
import numpy as np
import matplotlib.pyplot as plt
import glob as g
import os
gdal.UseExceptions()
Cargar una escena Landsat, unir todas las bandas en un arreglo y crear un nuevo raster que tenga todas las bandas en un sólo archivo.
In [2]:
landsat_tiff = "/Users/RaulSierra/Ecoinformatica/Divulgacion/Talleres/LUCC_ITAM/datos/ejercicio2.tif"
raster = gdal.Open(landsat_tiff)
print raster.GetGeoTransform()
ls_data = np.zeros((2391, 2429, raster.RasterCount))
for i in range(1, raster.RasterCount):
#Para ver la imagen
#Leemos cada banda al arreglo
ls_data[:,:,i-1] = raster.GetRasterBand(i).ReadAsArray()
Podemos hacer los mismo que hacemos en QGIS para visualizar (normalizar estadísticamente las bandas y aplastar todos lo valores a más de N desviaciones estándar a que estén exactamente a N*std. En este ejemplo, acotamos a 2 desviaciones estándar.
In [3]:
num_bands = raster.RasterCount - 1
norm_bands = ls_data[:,:,0:num_bands].copy()
print norm_bands.shape
means = norm_bands.reshape(-1, num_bands).mean(0)
sds = norm_bands.reshape(-1, num_bands).std(0)
print means
print sds
norm_bands -= means
norm_bands /= 2*sds
norm_bands[norm_bands > 1] = 1
norm_bands[norm_bands < -1] = -1
print norm_bands.reshape(-1, num_bands).mean(0)
print norm_bands.reshape(-1, num_bands).std(0)
norm_bands -= norm_bands.reshape(-1, num_bands).min(0)
norm_bands /= norm_bands.reshape(-1, num_bands).max(0)
print "Ahora los valores estan entre 0 y 1 para imshow"
print norm_bands.reshape(-1, num_bands).min(0)
print norm_bands.reshape(-1, num_bands).max(0)
Y visualizamos con color "real"
In [4]:
plt.axis('off')
plt.imshow(norm_bands[:,:,(2, 1, 0)])
Out[4]:
O con color falso para ver vegetación en rojo:
In [5]:
plt.axis('off')
plt.imshow(norm_bands[:,:,(3, 4, 2)])
Out[5]:
Vamos a convertir el atributo que corresponde al nivel 1 de la clasificación de coberturas de suelo de cadena de caracteres a numérico, con el objetivo de generar un raster que podamos procesar con las imágenes landsat.
Primero abrimos y leemos el shapefile fuente.
In [18]:
import fiona
shape_dir = "../datos/interpretados/CHIPS_AGUASCALIENTES_LANDSAT_2000/"
with fiona.open(os.path.join(shape_dir, "CHIPS_AGUASCALIENTES_LANDSAT_2000.shp")) 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'])
Luego creamos un nuevo shapefile y mapeamos los nombres en el atributo level_1 a valores numéricos
In [7]:
#Queremos que el nivel 1 sea numérico en vez de una cadena de caracteres
src_schema['properties']
with fiona.open(os.path.join(shape_dir, "CHIPS_AGUASCALIENTES_LANDSAT_2000_L1.shp"),
'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
Out[7]:
Para rasterizar el shapefile usamos ogr, que viene junto con gdal en la biblioteca osgeo.
Primero abrimos el shapefile que creamos en el paso anterior.
In [8]:
from osgeo import ogr
interpretados = ogr.Open("../datos/interpretados/CHIPS_AGUASCALIENTES_LANDSAT_2000/CHIPS_AGUASCALIENTES_LANDSAT_2000_L1.shp")
# No confundir features de shapefiles con características en machine learning
layer = interpretados.GetLayer()
featureCount = layer.GetFeatureCount()
print featureCount
lyrDefn = layer.GetLayerDefn()
# Print campos del shapefile
for i in range(lyrDefn.GetFieldCount()):
print lyrDefn.GetFieldDefn(i).GetName()
Tenemos que definir cómo convertir de polígonos a pixeles. Esto se hace definiendo el tamaño que ocupa un pixel de la región que define el shapefile.
In [13]:
# Define pixel_size and NoData value of new raster
pixel_size = 30
NoData_value = 0
# Open the data source and read in the extent
source_srs = layer.GetSpatialRef()
x_min, x_max, y_min, y_max = layer.GetExtent()
print [x_min, x_max, y_min, y_max]
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
print [x_res, y_res]
target_ds = gdal.GetDriverByName('GTiff').Create('../datos/interpretados/CHIPS_AGUASCALIENTES_LANDSAT_2000/CHIPS_AGUASCALIENTES_LANDSAT_2000_L1.tiff', x_res, y_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)
gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=level_1"])
if source_srs:
print "Aplicar proyeccion"
# Make the target raster have the same projection as the source
target_ds.SetProjection(source_srs.ExportToWkt())
In [14]:
import seaborn as sns
cpal = sns.choose_colorbrewer_palette("qualitative", as_cmap=False)
In [15]:
import matplotlib.cm as cm
clases_im = band.ReadAsArray()
cpal[8] = (1, 1, 1)
snspal = cm.Accent.from_list("Accent", cpal)
clases_im[clases_im == 0] = 8 #para que ponga los nodata blanco
#cerramos el raster
target_ds = None
plt.axis('off')
plt.imshow(clases_im, cmap = snspal)
Out[15]:
Con el siguiente código vamos a crear un raster que le asigne IDs a los pixeles de nuestra imagen de tal forma que los pixeles con el mismo ID pertenezcan a la misma área y tengan la misma clase.
Para esto vamos a crear un shapefile a partir de la capa de etiquetas (banda 10) de nuestro raster (ejercicio2.tif). Utilizamos las herramientas de línea de comandos de gdal.
In [ ]:
In [ ]:
src_file = "../datos/ejercicio2.tif"
dst_file = "../datos/ejercicio2_polIDS"
#este comando lo pueden correr directamente en la consola
!gdal_polygonize.py $src_file -f "ESRI Shapefile" -b 10 $dst_file ejercicio2_polIDS Etiqueta
El shapefile que obtenemos está compuesto por polígonos que corresponden a conjuntos de pixeles conectados en el raster que comparten la misma etiqueta, en este nuevo shapefile los polígonos no tienen ID así que les asigaremos a cada uno un ID único en un nuevo campo "OID".
No encontré cómo hacer lo anterior con las herramientas de línea de comando así que lo haré en python.
In [ ]:
ej_shp = ogr.Open("../datos/ejercicio2_polIDS/ejercicio2_polIDS.shp", 1)
# No confundir features de shapefiles con características en machine learning
layer = ej_shp.GetLayer()
featureCount = layer.GetFeatureCount()
print featureCount
idField = ogr.FieldDefn("OID", ogr.OFTInteger)
layer.CreateField(idField)
for i in range(featureCount):
f = layer.GetFeature(i)
f.SetField("OID", i)
layer.SetFeature(f)
ej_shp = None
El código anterior actualizó el shapefile con el nuevo campo 'OID'. Ahora rasterizamos el shape escribiendo los valores OID en los pixeles.
In [ ]:
from osgeo import gdal
src_file = "../datos/ejercicio2_polIDS/ejercicio2_polIDS.shp"
dst_file = "../datos/ejercicio2_polIDS.tif"
raster = gdal.Open("../datos/ejercicio2.tif")
width = raster.RasterXSize
height = raster.RasterYSize
!gdal_rasterize -a OID -a_nodata 0 -ts $width $height -l ejercicio2_polIDS $src_file $dst_file
In [9]:
cnt = gdal.GetDriverCount()
formatsList = [] # Empty List
for i in range(cnt):
driver = gdal.GetDriver(i)
driverName = driver.ShortName
driverLong = driver.LongName
if not driverName in formatsList:
formatsList.append((driverName, driverLong))
formatsList.sort() # Sorting the messy list of ogr drivers
for i in formatsList:
print i