Преобразовать файл сетки ASCII в GeoTIFF, используя Python?

11

У меня есть файл формата сетки ASCII. Например:

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

Как я могу преобразовать его в TIFF или любой другой растр, используя Python?

voimak
источник
Программное обеспечение ГИС может конвертировать аси в геотифф. Кодирование не требуется. Я использую QGIS. Это бесплатно.
Сол Шеард

Ответы:

13

Версия с псевдокодом:

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.

Образец, который поможет вам в этом - отсюда :

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)
Джей Лаура
источник
а значения 14,97 и -34,54 являются координатами верхнего левого угла в координатах WGS84?
Слава
9

Альтернатива с использованием gdal_translate :

gdal_translate -of "GTiff" fname.asc outname.tif
bananafish
источник
7

Создание копии может быть проще, поскольку ваш файл представляет собой AAIGrid, а GTiff поддерживает 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

Любой драйвер, который поддерживает CreateCopy, может использовать это.


источник
Если вам не нужно использовать Python, Bananafish определенно подходит.
удивительно, спасибо! Мой входной файл .asc поставляется без CRS. Есть ли способ указать этот CRS (GCS WGS 84) перед записью растра?
RutgerH
используйте SetProjection и строку. Вы можете получить строку из OSR. Смотрите ответ редактировать.