Converter um arquivo de grade ASCII para GeoTIFF usando Python?

11

Eu tenho um arquivo de formato raster de grade ASCII. Por exemplo:

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

Como posso convertê-lo para TIFF ou qualquer outra varredura usando Python?

voimak
fonte
O software GIS pode converter asci em geotiff. Nenhuma codificação necessária. Eu uso o QGIS. É grátis.
Saul Sheard

Respostas:

13

A versão do pseudo-có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.

Uma amostra que irá ajudá-lo - a partir daqui :

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
fonte
e os valores 14,97 e -34,54 são as coordenadas do canto superior esquerdo nos cooridanatos WGS84?
Slava
9

Alternativa usando gdal_translate :

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

Pode ser mais fácil fazer uma cópia de criação, pois seu arquivo é um AAIGrid e o GTiff suporta 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

Qualquer driver que ofereça suporte ao CreateCopy pode usar isso.


fonte
Se você não precisa usar python, o bananafish está definitivamente certo.
incrível, obrigado! Meu arquivo .asc de entrada vem sem um CRS. Existe uma maneira de especificar esse CRS (GCS WGS 84) antes de escrever a varredura?
perfil completo de Rutger
use SetProjection e uma string. Você pode obter a string do osr. Veja a edição de resposta.