[gdal-dev] saving numpy array as ascii raster

[Bryan Keith]

> Hi,
> Thanks for your quick reply.
> When I edit my code function to be:
> def WriteRaster (dst_filename, raster):
>
>         format = "AAIGrid"
>         driver = gdal.GetDriverByName( format )
>         dst_ds = driver.Create( dst_filename, 71, 73,\
>               1,gdal.GDT_Float32,options=["COMPRESS=PACKBITS","TFW=YES"] )
>
> Python complains:
>
> ERROR 6: GDALDriver::Create() ... no create method implemented for this
> format.
>

Oz,

AAIGrid cannot use Create, only CreateCopy.  That's the difference between
the + in the gdalinfo --formats list.

I tried this code, and I think it does what you're looking for:

#! /usr/bin/python
from osgeo import gdal
from osgeo.gdalconst import *
dataset = gdal.Open('raster1.asc', GA_ReadOnly)
dataset1 = gdal.Open('raster2.asc', GA_ReadOnly)
print dataset, dataset1
from osgeo import gdal_array
from osgeo import osr
a = gdal_array.DatasetReadAsArray(dataset)
b = gdal_array.DatasetReadAsArray(dataset1)
print a
print b
def WriteRaster (dst_filename, raster):

        format = "MEM"
        driver = gdal.GetDriverByName( format )
        dst_ds = driver.Create( dst_filename, len(raster[0]), len(raster),\
               1,gdal.GDT_Float32)
        dst_ds.SetGeoTransform( [-19.5, 1.0, 0.0, 37.5, 0.0, -1.0] )
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326) #WGS84 lat long.
        dst_ds.SetProjection( srs.ExportToWkt() )
        dst_ds.GetRasterBand(1).WriteArray( raster )


        format = 'AAIGrid'
        driver = gdal.GetDriverByName(format)
        dst_ds_new = driver.CreateCopy(dst_filename, dst_ds)

        dst_ds = None
c= a+b
print c
WriteRaster ('raster3.asc', c)

Bryan