Convertir un fichier de grille ASCII en GeoTIFF en utilisant Python?

11

J'ai un fichier au format raster de grille ASCII. Par exemple:

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

Comment puis-je le convertir en TIFF ou tout autre raster en utilisant Python?

voimak
la source
Le logiciel SIG peut convertir asci en géotiff. Aucun codage nécessaire. J'utilise QGIS. C'est gratuit.
Saul Sheard

Réponses:

13

La version du pseudo code:

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.

Un échantillon qui vous aidera - d'ici :

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
la source
et les valeurs 14,97 et -34,54 sont les coordonnées du coin supérieur gauche dans les cooridanates WGS84?
Slava
9

Alternative utilisant gdal_translate :

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

Il peut être plus facile de créer une copie, car votre fichier est un AAIGrid et GTiff prend en charge 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

Tout pilote prenant en charge CreateCopy peut l'utiliser.


la source
Si vous n'avez pas besoin d'utiliser python, le bananafish a définitivement raison.
incroyable, merci! Mon fichier .asc d'entrée est livré sans CRS. Existe-t-il un moyen de spécifier ce CRS (GCS WGS 84) avant d'écrire le raster?
RutgerH
utilisez SetProjection et une chaîne. Vous pouvez obtenir la chaîne d'osr. Voir réponse éditer.