Konwertuj plik siatki ASCII na GeoTIFF za pomocą Pythona?

11

Mam plik formatu rastrowego ASCII. Na przykład:

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

Jak mogę przekonwertować go na TIFF lub inny raster za pomocą Pythona?

voimak
źródło
Oprogramowanie GIS może konwertować ASCI na geotiff. Kodowanie nie jest konieczne. Używam QGIS. Jest wolne.
Saul Sheard

Odpowiedzi:

13

Wersja pseudokodu:

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.

Próbka, która pomoże ci dalej - stąd :

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
źródło
a wartości 14.97 i -34.54 są współrzędnymi lewego górnego rogu w współrzędnych WGS84?
Slava
9

Alternatywne użycie gdal_translate :

gdal_translate -of "GTiff" fname.asc outname.tif
bananowiec
źródło
7

Utworzenie kopii może być łatwiejsze, ponieważ plik ma format AAIGrid, a GTiff obsługuje funkcję 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

Może to wykorzystać dowolny sterownik obsługujący CreateCopy.


źródło
Jeśli nie musisz używać pytona, bananowiec ma zdecydowanie rację.
niesamowite, dzięki! Mój wejściowy plik .asc nie ma CRS. Czy istnieje sposób na określenie tego CRS (GCS WGS 84) przed zapisaniem rastra?
RutgerH
użyj SetProjection i łańcucha. Możesz pobrać ciąg z osr. Zobacz edycję odpowiedzi.