[Gvsig_usuarios] Transformação radiométrica de imagens Landsat 8 (16 bits para 8 bits)

[Óscar Martínez]

Buenas,

He realizado un pequeño script que realizaría la transformación de 16 
bits a 8 bits usando solo la librería gdal. Es como usar gdal_translate 
pero en forma de código.

Aquí explica como cargar el script:

https://www.youtube.com/watch?v=7c_6KetDOAM

Y el script lo adjunto como texto al correo.

Ten una Vista abierta antes de ejecutarlo. Al ejecutar el script solo te 
pedirá que selecciones la ruta del fichero, y al acabar la imagen 
resultado se cargará automáticamente en la Vista.

Por las pruebas que he realizado creo que funciona correctamente la 
coversión entre 16 a 8 bits. Si puedes probarlo y comentarme resultados 
y si te deja calcular el ndvi estaría genial.

Si tienes alguna sugerencia podría mejorar el script (como dar opción de 
seleccionar capas raster ya cargadas en una vista), si todo es correcto 
podría hacer que se pueda incorporar como geoproceso. Pero primero me 
gustaría saber si os es de utilidad.

Puedes cargar imagenes de 16 bits de Landsat y sentinel si no me 
equivoco sin problemas (es algo que me ha comentado Gilberto), pero para 
el cálculo ndvi no se exacto que necesitas.

Prueba y me comentas :)

Un saludo


El 28/10/16 a las 14:49, Eliazer Kosciuk escribió:
> Buenas, pessoal!
>
> Em primeiro lugar, desculpem-me por escrever em português, mas tenho 
> certeza que serei melhor entendido do que se eu tentar escrever no meu 
> "portunhol". Podem responder em espanhol...
>
> Estava tentando usar o gvSIG 2.3 para gerar um NDVI com imagens 
> Landsat 8, e o gvSIG simplesmente gera uma imagem sem informações. 
> Pesquisando e perguntando por aí, descobri que o gvSIG ainda não lida 
> com as imagens no formato 16 bits, que é o formato entregue pelas 
> imagens Landsat 8.
>
> Pesquisando sobre como efetuar a conversão de uma imagem 16 bits para 
> 8 bits, encontrei esse tutorial do Jorge Santos, para o QGIS: 
> http://www.processamentodigital.com.br/2014/12/12/qgis-2-6-transformacao-radiometrica-para-imagens-landsat-8-16-bit-para-8-bit/
>
> Considerando que o gvSIG também contém a biblioteca GDAL, há a 
> possibilidade de realizar esse procedimento no gvSIG? Ou existe algum 
> outro caminho para trabalhar com imagens Landsat8 no gvSIG, para 
> geração de NDVI?
>
> [ ]s,
>
> Eliazer Kosciuk
> http://geo.ideaplus.com.br
>
>
>
>
> _______________________________________________
> gvSIG_usuarios mailing list
> gvSIG_usuarios en listserv.gva.es
>
> Para ver histórico de mensajes, editar sus preferencias de usuario o darse de baja en esta lista, acuda a la siguiente dirección:
>
> https://listserv.gva.es/cgi-bin/mailman/listinfo/gvsig_usuarios

------------ próxima parte ------------
Se ha borrado un adjunto en formato HTML...
URL: <http://listserv.gva.es/pipermail/gvsig_usuarios/attachments/20161031/bf89097e/attachment.html>
------------ próxima parte ------------
 
# encoding: utf-8


from gvsig import *
from gvsig import uselib
uselib.use_plugin("org.gvsig.gdal.app.mainplugin")
from gvsig import geom
from org.gdal.gdal import gdal
import os
from org.gdal.gdalconst import gdalconst
from gvsig import commonsdialog

import threading

def main(*args):
        threading.Thread(target=gdal16to8, name="Gdal", args=tuple()).start()

def gdal16to8(*args):
    uselib.use_plugin("org.gvsig.gdal.app.mainplugin")
    from gvsig import geom
    from org.gdal.gdal import gdal
    import os
    from org.gdal.gdalconst import gdalconst
    import jarray
    
    #sourcefile = "/home/osc/temp/gdal/LC81980332016298LGN00_B6.TIF"
    sourcefile = commonsdialog.openFileDialog(title='', initialPath=None, root=None)[0]
    gdal.AllRegister()
    
    if not os.path.exists(sourcefile):
        print "- capa no encontrada"
        return
    import ntpath
    pathfile = ntpath.basename(sourcefile)
    namefile = os.path.splitext(pathfile)[0]

    #output_file = getTempFile("/home/osc/temp/gdal/r", ".TIF")
    output_file = getTempFile("16to8-"+str(namefile), ".TIF")
    print output_file
    
    ### GDAL: Read Raster
    dataset = gdal.Open(sourcefile)
    print "dataset: ", dataset, type(dataset)
     
    
    geotransform = dataset.GetGeoTransform()
    prj=dataset.GetProjection()

    #values
    
    band = dataset.GetRasterBand(1)
    #from java.lang import Double
    rmin = [None]#[Double(0.0)]
    rmax = [None]#[Double(0.0)]

    band.GetMinimum(rmin)
    band.GetMaximum(rmax)
  
    print "Pre min, max: ", rmin, rmax
    
    if rmin is [None] or rmax is [None]:
        (rmin,rmax) = band.ComputeRasterMinMax(1)
        
    print "ComputeRasterMinMax: ", rmin, rmax

    bmin = jarray.array([0.0], "d")
    bmax = jarray.array([0.0], "d")
    bmean = jarray.array([0.0], "d")
    bstv = jarray.array([0.0], "d")
    
    band.GetStatistics(True, True, bmin, bmax, bmean, bstv)

    print "\n== Statistics"
    print "Min: ", bmin[0]
    print "Max: ", bmax[0]
    print "Mean: ", bmean[0]
    print "Stv: ", bstv[0]
    
    rmin = int(bmin[0])
    rmax = int(bmax[0])
    print "rmin, rmax: ", rmin, rmax
    print type(rmin), type(rmax)
    
    
    xsize = band.getXSize()
    ysize = band.getYSize()
    print "Size x, y: ", xsize, ysize

    #xsize = 4000
    #ysize = 4000
    
    # Creating output
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(output_file, xsize, ysize, 1, gdalconst.GDT_Byte)

    # Writing partial
    for i in range(0, ysize):
        values = jarray.array([0]* xsize * 1, "i") # "d"
        band.ReadRaster(0, i, xsize, 1, values)
        #print values
        for a in range(0, len(values)):
            values[a] = values[a]*255/rmax

        dst_ds.GetRasterBand(1).WriteRaster(0, i, xsize, 1, values)

    dst_ds.SetGeoTransform( [ geotransform[0], geotransform[1], 0, geotransform[3], 0, geotransform[5] ] )

    # set projection of new raster
    dst_ds.SetProjection( prj )

    dst_ds.delete()
    
    r = loadRasterFile(output_file)
    r.setName(namefile+"_8bits")

    print "Done"