ASCII file with latitude, longitude and data to GeoTiff using Python

Question

I have a .asc file with three columns of data: longitude, latitude and a data value. Each latitude and longitude point represents the location of a centriod of a 1/8 degree unprojected grid. I want to convert this file to a raster (.tif) using gdal in Python. I then want to use projected datasets to subset this large grid. Any thoughts on how to begin here?

Answer 1

To convert an ASCII file with longitude, latitude and data value you may use a function like this:

from osgeo import gdal

def csv2tif(source, target):
    cvs = gdal.Open(source)
    if cvs is None:
        print 'ERROR: Unable to open %s' % source
        return

    geotiff = gdal.GetDriverByName("GTiff")
    if geotiff is None:
        print 'ERROR: GeoTIFF driver not available.'
        return

    options = []
    geotiff.CreateCopy(target, cvs, 0, options)

source = 'E:\\test.csv'
target = 'E:\\test.tif'

csv2tif(source, target)

The next function imports and reprojects the CSV and saves it as geotiff file. When you want the function to handle other formats as well, some more parameters are required. You also may change the algorithm (gdal.GRA_Bilinear) to calculate the cell values of the projected raster or change the its resolution (xsize, ysize).

from osgeo import gdal, osr

def csv2tif_projected(source, target, destEPSG, srcEPSG=4326):

    # open CSV source file
    cvs = gdal.Open(source)
    if cvs is None:
        print 'ERROR: Unable to open %s' % source
        return

    # get GeoTIFF driver
    geotiff = gdal.GetDriverByName("GTiff")
    if geotiff is None:
        print 'ERROR: GeoTIFF driver not available.'
        return

    # set source coordinate system of coordinates in CSV file
    src_crs = osr.SpatialReference()
    src_crs.ImportFromEPSG(srcEPSG)

    # set destination projection parameters
    dest_crs = osr.SpatialReference()
    dest_crs.ImportFromEPSG(destEPSG)

    # set coordinate transformation
    tx = osr.CoordinateTransformation(src_crs, dest_crs)

    # get raster dimension related parameters of source dataset
    xo, xs, xr, yo, yr, ys = cvs.GetGeoTransform()
    xsize = cvs.RasterXSize
    ysize = cvs.RasterYSize

    # convert corner coordinates from old to new coordinate system
    (ulx, uly, ulz) = tx.TransformPoint(xo, yo)
    (lrx, lry, lrz) = tx.TransformPoint(xo + xs * xsize + xr * ysize,\
                                        yo + yr * xsize + ys * ysize)

    # create blank in-memory raster file with same dimension as CSV raster
    mem = gdal.GetDriverByName('MEM')
    dest_ds = mem.Create('', xsize, ysize, 1, gdal.GDT_Float32)

    # get new transformation
    dest_geo = (ulx, (lrx - ulx) / xsize, xr,\
                uly, yr, (lry - uly) / ysize)

    # set the geotransformation
    dest_ds.SetGeoTransform(dest_geo)
    dest_ds.SetProjection(dest_crs.ExportToWkt())

    # project the source raster to destination coordinate system
    gdal.ReprojectImage(cvs, dest_ds, \
                        src_crs.ExportToWkt(), dest_crs.ExportToWkt(),\
                        gdal.GRA_Bilinear, 0.0, 0.0)

    # save projected in-memory raster to disk
    geotiff.CreateCopy(target, dest_ds, 0 )

For clipping raster with shapefile see the explanation and code in Python GDAL/OGR Cookbook Clip a GeoTiff with Shapefile as well as the discussion Clipping raster with vector layer using GDAL.

Answer 2

I ended up building a function that uses a reference image with the correct projection and extent to transform the XYZ ascii file into a GeoTiff. My ascii file had no header, so gdal_translate didn't work. Here's the function:

def ascii_to_tiff(infile, outfile, refIm):
    """
    Transform an XYZ ascii file without a header to a projected GeoTiff

    :param infile (str): path to infile ascii location
    :param outfile (str): path to final GTiff
    :param refIm (str): path to a reference image made from the lat lon pair centriods

    """

    im = gdal.Open(refIm)
    ima = gdal.Open(refIm).ReadAsArray()
    row = ima.shape[0];
    col = ima.shape[1]

    indata = np.genfromtxt(infile, delimiter=",", skip_header=True, dtype=None)
    lon = indata[:, 0]  # x
    lat = indata[:, 1]  # y
    data = indata[:, 2]

    # create grid
    xmin, xmax, ymin, ymax = [min(lon), max(lon), min(lat), max(lat)]
    xi = np.linspace(xmin, xmax, col)
    yi = np.linspace(ymin, ymax, row)
    xi, yi = np.meshgrid(xi, yi)

    # linear interpolation
    zi = ml.griddata(lon, lat, data, xi, yi, interp='linear')
    final_array = np.asarray(np.rot90(np.transpose(zi)))

    # projection
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(outfile, col, row, 1, gdal.GDT_Float32)
    dst_ds.GetRasterBand(1).WriteArray(final_array)
    prj = im.GetProjection()
    dst_ds.SetProjection(prj)

    gt = im.GetGeoTransform()
    dst_ds.SetGeoTransform(gt)
    dst_ds = None

    final_tif = gdal.Open(outfile, GA_ReadOnly).ReadAsArray()
    return final_tif

Answer 3

Using the core GDAL utilities, you could easily convert your ASCII file to GTiff using gdal_translate as follows:

gdal_translate -of 'GTiff' input.asc output.tiff

Answer 4

Wrap up another answer with Python using os.system or subprocess. Although subprocess should be preferred for various reasons, see here

import os, subprocess

os.system("gdal_translate -of 'GTiff' input.asc output.tiff")

# or
subprocess.check_ouput(["gdal_translate", "-of GTiff", "input.asc", "output.tif"])