¿Convertir un archivo de cuadrícula ASCII a GeoTIFF usando Python?

11

Tengo un archivo de formato de trama de cuadrícula ASCII. Por ejemplo:

ncols 480
nrows 450
xllcorner 378923
yllcorner 4072345
cellsize 30
nodata_value -32768
43 2 45 7 3 56 2 5 23 65 34 6 32 54 57 34 2 2 54 6 
35 45 65 34 2 6 78 4 2 6 89 3 2 7 45 23 5 8 4 1 62 ...

¿Cómo puedo convertirlo a TIFF o cualquier otro ráster usando Python?

python gdal convert voimak
fuente
El software SIG puede convertir asci a geotiff. No se necesita codificación. Yo uso QGIS. Es gratis.
Saul Sheard

Respuestas:

13

La versión del pseudocódigo:

import gdal
import numpy

create the gdal output file as geotiff
set the no data value
set the geotransform 

numpy.genfromtxt('your file', numpy.int8) #looks like int from you example
reshape your array to the shape you need

write out the array.

Una muestra que lo ayudará a lo largo, desde aquí :

if __name__ == '__main__':
    # Import libs
    import numpy, os
    from osgeo import osr, gdal

    # Set file vars
    output_file = "out.tif"

    # Create gtif
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(output_file, 174, 115, 1, gdal.GDT_Byte )
    raster = numpy.zeros( (174, 115) )

    # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
    dst_ds.SetGeoTransform( [ 14.97, 0.11, 0, -34.54, 0, 0.11 ] )

    # set the reference info 
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS("WGS84")
    dst_ds.SetProjection( srs.ExportToWkt() )

    # write the band
    dst_ds.GetRasterBand(1).WriteArray(raster)
Jay Laura
fuente
y los valores 14.97 y -34.54 son las coordenadas de la esquina superior izquierda en las coordenadas WGS84?
Slava
9

Alternativa usando gdal_translate :

gdal_translate -of "GTiff" fname.asc outname.tif
bananafish
fuente
7

Puede ser más fácil hacer una copia de creación, ya que su archivo es AAIGrid y GTiff admite CreateCopy ():

from osgeo import gdal, osr
drv = gdal.GetDriverByName('GTiff')
ds_in = gdal.Open('in.asc')
ds_out = drv.CreateCopy('out.tif', ds_in)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds_out.SetProjection(srs.ExportToWkt())
ds_in = None
ds_out = None

Cualquier controlador que admita CreateCopy puede usar esto.


fuente
Si no necesita usar python, bananafish definitivamente tiene razón.
increíble, gracias! Mi archivo de entrada .asc viene sin un CRS. ¿Hay alguna manera de especificar este CRS (GCS WGS 84) antes de escribir el ráster?
RutgerH
use SetProjection y una cadena. Puede obtener la cadena de osr. Ver respuesta editar.