"""Interpolate csv output results to raster format
Jonathan Griffin
Geoscience Australia June 2017
"""

import os, sys
from subprocess import call
from glob import glob
import numpy as np
from scipy.interpolate import griddata
import ogr, gdal, osr

def grid_results(infile, resolution = 0.01, clip_shp = None, 
                 overwrite=True, contour=False):
    """Read csv file, interpolate and convert to raster
    """
    outfile = infile.rstrip('().csv') + '_gridded.tif'
   # if not overwrite:
    if os.path.isfile(outfile):
        if not overwrite:
            print('Not creating file %s as already exists' % outfile)
            print('To re-create file (e.g if inputs changed) set overwrite=True)')
            return
        else:
            try:
                os.remove(outfile)
                os.remove((outfile.rstrip('.tif') + '_clip.tif'))
            except:
                pass
    data = np.genfromtxt(infile, delimiter=',')
    max_lon = max(data[:,0])
    min_lon = min(data[:,0])
    max_lat = max(data[:,1])
    min_lat = min(data[:,1])
    #print max_lon, min_lon, max_lat, min_lat
    xi = np.arange(min_lon, max_lon, resolution)
    yi = np.arange(min_lat, max_lat, resolution)
    XI,YI = np.meshgrid(xi,yi)
    xsize = len(xi)
    ysize = len(yi)

    print('Interpolating results')
    gridded_results = griddata((data[:,0],data[:,1]),data[:,2],(XI,YI),method='linear')
    #print gridded_results
    #outfile = infile.rstrip('().csv') + '_gridded.tif'
    print('Writing gridded data to %s' % outfile)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Float32)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    gt = [(min_lon - (resolution/2)), resolution, 0, 
          (min_lat - (resolution/2)), 0, resolution]
    ds.SetGeoTransform(gt)
    outband=ds.GetRasterBand(1)
    outband.SetStatistics(np.min(gridded_results), np.max(gridded_results), np.average(gridded_results), np.std(gridded_results))
    outband.WriteArray(gridded_results)
    # Need to close output dataset before we can do clipping
    ds = None
    # now clip by shapefile
    if clip_shp is not None:
        clipfile = outfile.rstrip('.tif') + '_clip.tif'
        cmd = ['gdalwarp',
               '-cutline',
               clip_shp,
               '-crop_to_cutline',
               '-dstalpha',
               outfile,
               clipfile]
        print(cmd)
        call(cmd, shell=False)
    if contour is True:
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (outfile,  outfile.rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (clipfile, clipfile.rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
   # outfile = infile.rstrip('.csv') + '_gridded.csv'
   # np.savetxt(outfile, gridded_results)

filelist = glob('outputs/*mmi*1780*sq*Chiou*.csv')
#infile = 'outputs/scenario_gmf_loc_mmi_1699megathrust_AtkinsonBoore2003SInter().csv'
clipping_shapefile = '/short/n74/jdg547/eq_hazmap_tests/data/IDN_Admin/IDN_adm0.shp'
for infile in filelist:
    grid_results(infile, resolution=0.01, clip_shp = clipping_shapefile,
                 overwrite=True, contour=True)

42f62fcf349d5fc05b05f3b193b75ac956b5c024
<mask token>

<mask token>


def grid_results(infile, resolution=0.01, clip_shp=None, overwrite=True,
    contour=False):
    """Read csv file, interpolate and convert to raster
    """
    outfile = infile.rstrip('().csv') + '_gridded.tif'
    if os.path.isfile(outfile):
        if not overwrite:
            print('Not creating file %s as already exists' % outfile)
            print(
                'To re-create file (e.g if inputs changed) set overwrite=True)'
                )
            return
        else:
            try:
                os.remove(outfile)
                os.remove(outfile.rstrip('.tif') + '_clip.tif')
            except:
                pass
    data = np.genfromtxt(infile, delimiter=',')
    max_lon = max(data[:, 0])
    min_lon = min(data[:, 0])
    max_lat = max(data[:, 1])
    min_lat = min(data[:, 1])
    xi = np.arange(min_lon, max_lon, resolution)
    yi = np.arange(min_lat, max_lat, resolution)
    XI, YI = np.meshgrid(xi, yi)
    xsize = len(xi)
    ysize = len(yi)
    print('Interpolating results')
    gridded_results = griddata((data[:, 0], data[:, 1]), data[:, 2], (XI,
        YI), method='linear')
    print('Writing gridded data to %s' % outfile)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Float32)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    gt = [min_lon - resolution / 2, resolution, 0, min_lat - resolution / 2,
        0, resolution]
    ds.SetGeoTransform(gt)
    outband = ds.GetRasterBand(1)
    outband.SetStatistics(np.min(gridded_results), np.max(gridded_results),
        np.average(gridded_results), np.std(gridded_results))
    outband.WriteArray(gridded_results)
    ds = None
    if clip_shp is not None:
        clipfile = outfile.rstrip('.tif') + '_clip.tif'
        cmd = ['gdalwarp', '-cutline', clip_shp, '-crop_to_cutline',
            '-dstalpha', outfile, clipfile]
        print(cmd)
        call(cmd, shell=False)
    if contour is True:
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (outfile, outfile.
            rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (clipfile, clipfile.
            rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)


<mask token>
for infile in filelist:
    grid_results(infile, resolution=0.01, clip_shp=clipping_shapefile,
        overwrite=True, contour=True)

<mask token>


def grid_results(infile, resolution=0.01, clip_shp=None, overwrite=True,
    contour=False):
    """Read csv file, interpolate and convert to raster
    """
    outfile = infile.rstrip('().csv') + '_gridded.tif'
    if os.path.isfile(outfile):
        if not overwrite:
            print('Not creating file %s as already exists' % outfile)
            print(
                'To re-create file (e.g if inputs changed) set overwrite=True)'
                )
            return
        else:
            try:
                os.remove(outfile)
                os.remove(outfile.rstrip('.tif') + '_clip.tif')
            except:
                pass
    data = np.genfromtxt(infile, delimiter=',')
    max_lon = max(data[:, 0])
    min_lon = min(data[:, 0])
    max_lat = max(data[:, 1])
    min_lat = min(data[:, 1])
    xi = np.arange(min_lon, max_lon, resolution)
    yi = np.arange(min_lat, max_lat, resolution)
    XI, YI = np.meshgrid(xi, yi)
    xsize = len(xi)
    ysize = len(yi)
    print('Interpolating results')
    gridded_results = griddata((data[:, 0], data[:, 1]), data[:, 2], (XI,
        YI), method='linear')
    print('Writing gridded data to %s' % outfile)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Float32)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    gt = [min_lon - resolution / 2, resolution, 0, min_lat - resolution / 2,
        0, resolution]
    ds.SetGeoTransform(gt)
    outband = ds.GetRasterBand(1)
    outband.SetStatistics(np.min(gridded_results), np.max(gridded_results),
        np.average(gridded_results), np.std(gridded_results))
    outband.WriteArray(gridded_results)
    ds = None
    if clip_shp is not None:
        clipfile = outfile.rstrip('.tif') + '_clip.tif'
        cmd = ['gdalwarp', '-cutline', clip_shp, '-crop_to_cutline',
            '-dstalpha', outfile, clipfile]
        print(cmd)
        call(cmd, shell=False)
    if contour is True:
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (outfile, outfile.
            rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (clipfile, clipfile.
            rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)


filelist = glob('outputs/*mmi*1780*sq*Chiou*.csv')
clipping_shapefile = (
    '/short/n74/jdg547/eq_hazmap_tests/data/IDN_Admin/IDN_adm0.shp')
for infile in filelist:
    grid_results(infile, resolution=0.01, clip_shp=clipping_shapefile,
        overwrite=True, contour=True)

<mask token>
import os, sys
from subprocess import call
from glob import glob
import numpy as np
from scipy.interpolate import griddata
import ogr, gdal, osr


def grid_results(infile, resolution=0.01, clip_shp=None, overwrite=True,
    contour=False):
    """Read csv file, interpolate and convert to raster
    """
    outfile = infile.rstrip('().csv') + '_gridded.tif'
    if os.path.isfile(outfile):
        if not overwrite:
            print('Not creating file %s as already exists' % outfile)
            print(
                'To re-create file (e.g if inputs changed) set overwrite=True)'
                )
            return
        else:
            try:
                os.remove(outfile)
                os.remove(outfile.rstrip('.tif') + '_clip.tif')
            except:
                pass
    data = np.genfromtxt(infile, delimiter=',')
    max_lon = max(data[:, 0])
    min_lon = min(data[:, 0])
    max_lat = max(data[:, 1])
    min_lat = min(data[:, 1])
    xi = np.arange(min_lon, max_lon, resolution)
    yi = np.arange(min_lat, max_lat, resolution)
    XI, YI = np.meshgrid(xi, yi)
    xsize = len(xi)
    ysize = len(yi)
    print('Interpolating results')
    gridded_results = griddata((data[:, 0], data[:, 1]), data[:, 2], (XI,
        YI), method='linear')
    print('Writing gridded data to %s' % outfile)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Float32)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    gt = [min_lon - resolution / 2, resolution, 0, min_lat - resolution / 2,
        0, resolution]
    ds.SetGeoTransform(gt)
    outband = ds.GetRasterBand(1)
    outband.SetStatistics(np.min(gridded_results), np.max(gridded_results),
        np.average(gridded_results), np.std(gridded_results))
    outband.WriteArray(gridded_results)
    ds = None
    if clip_shp is not None:
        clipfile = outfile.rstrip('.tif') + '_clip.tif'
        cmd = ['gdalwarp', '-cutline', clip_shp, '-crop_to_cutline',
            '-dstalpha', outfile, clipfile]
        print(cmd)
        call(cmd, shell=False)
    if contour is True:
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (outfile, outfile.
            rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (clipfile, clipfile.
            rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)


filelist = glob('outputs/*mmi*1780*sq*Chiou*.csv')
clipping_shapefile = (
    '/short/n74/jdg547/eq_hazmap_tests/data/IDN_Admin/IDN_adm0.shp')
for infile in filelist:
    grid_results(infile, resolution=0.01, clip_shp=clipping_shapefile,
        overwrite=True, contour=True)

"""Interpolate csv output results to raster format
Jonathan Griffin
Geoscience Australia June 2017
"""

import os, sys
from subprocess import call
from glob import glob
import numpy as np
from scipy.interpolate import griddata
import ogr, gdal, osr

def grid_results(infile, resolution = 0.01, clip_shp = None, 
                 overwrite=True, contour=False):
    """Read csv file, interpolate and convert to raster
    """
    outfile = infile.rstrip('().csv') + '_gridded.tif'
   # if not overwrite:
    if os.path.isfile(outfile):
        if not overwrite:
            print('Not creating file %s as already exists' % outfile)
            print('To re-create file (e.g if inputs changed) set overwrite=True)')
            return
        else:
            try:
                os.remove(outfile)
                os.remove((outfile.rstrip('.tif') + '_clip.tif'))
            except:
                pass
    data = np.genfromtxt(infile, delimiter=',')
    max_lon = max(data[:,0])
    min_lon = min(data[:,0])
    max_lat = max(data[:,1])
    min_lat = min(data[:,1])
    #print max_lon, min_lon, max_lat, min_lat
    xi = np.arange(min_lon, max_lon, resolution)
    yi = np.arange(min_lat, max_lat, resolution)
    XI,YI = np.meshgrid(xi,yi)
    xsize = len(xi)
    ysize = len(yi)

    print('Interpolating results')
    gridded_results = griddata((data[:,0],data[:,1]),data[:,2],(XI,YI),method='linear')
    #print gridded_results
    #outfile = infile.rstrip('().csv') + '_gridded.tif'
    print('Writing gridded data to %s' % outfile)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Float32)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    gt = [(min_lon - (resolution/2)), resolution, 0, 
          (min_lat - (resolution/2)), 0, resolution]
    ds.SetGeoTransform(gt)
    outband=ds.GetRasterBand(1)
    outband.SetStatistics(np.min(gridded_results), np.max(gridded_results), np.average(gridded_results), np.std(gridded_results))
    outband.WriteArray(gridded_results)
    # Need to close output dataset before we can do clipping
    ds = None
    # now clip by shapefile
    if clip_shp is not None:
        clipfile = outfile.rstrip('.tif') + '_clip.tif'
        cmd = ['gdalwarp',
               '-cutline',
               clip_shp,
               '-crop_to_cutline',
               '-dstalpha',
               outfile,
               clipfile]
        print(cmd)
        call(cmd, shell=False)
    if contour is True:
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (outfile,  outfile.rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
        cmd = 'gdal_contour -i 1 -off 0.5 %s %s.shp' % (clipfile, clipfile.rstrip('.tif'))
        print(cmd)
        call(cmd, shell=True)
   # outfile = infile.rstrip('.csv') + '_gridded.csv'
   # np.savetxt(outfile, gridded_results)

filelist = glob('outputs/*mmi*1780*sq*Chiou*.csv')
#infile = 'outputs/scenario_gmf_loc_mmi_1699megathrust_AtkinsonBoore2003SInter().csv'
clipping_shapefile = '/short/n74/jdg547/eq_hazmap_tests/data/IDN_Admin/IDN_adm0.shp'
for infile in filelist:
    grid_results(infile, resolution=0.01, clip_shp = clipping_shapefile,
                 overwrite=True, contour=True)
