def test_gdal_contour_2():
    if test_cli_utilities.get_gdal_contour_path() is None:
        return 'skip'

    try:
        os.remove('tmp/contour.shp')
    except:
        pass
    try:
        os.remove('tmp/contour.dbf')
    except:
        pass
    try:
        os.remove('tmp/contour.shx')
    except:
        pass

    # put -3d just after -fl to test #2793
    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -fl 10 20 25 -3d tmp/gdal_contour.tif tmp/contour.shp')

    size = 160
    precision = 1. / size

    ds = ogr.Open('tmp/contour.shp')

    expected_envelopes = [ [ 1.25, 1.75, 49.25, 49.75 ],
                           [ 1.25+0.125, 1.75-0.125, 49.25+0.125, 49.75-0.125 ],
                           [ 1.25+0.125+0.0625, 1.75-0.125-0.0625, 49.25+0.125+0.0625, 49.75-0.125-0.0625 ] ]
    expected_height = [ 10, 20, 25 ]

    lyr = ds.ExecuteSQL("select * from contour order by elev asc")

    if lyr.GetFeatureCount() != len(expected_envelopes):
        print('Got %d features. Expected %d' % (lyr.GetFeatureCount(), len(expected_envelopes)))
        return 'fail'

    i = 0
    feat = lyr.GetNextFeature()
    while feat is not None:
        if feat.GetGeometryRef().GetZ(0) != expected_height[i]:
            print('Got %f as z. Expected %f' % (feat.GetGeometryRef().GetZ(0), expected_height[i]))
            return 'fail'
        envelope = feat.GetGeometryRef().GetEnvelope()
        if feat.GetField('elev') != expected_height[i]:
            print('Got %f. Expected %f' % (feat.GetField('elev'), expected_height[i]))
            return 'fail'
        for j in range(4):
            if abs(expected_envelopes[i][j] - envelope[j]) > precision/2*1.001:
                print('i=%d, wkt=%s' % (i, feat.GetGeometryRef().ExportToWkt()))
                print(feat.GetGeometryRef().GetEnvelope())
                print(expected_envelopes[i])
                print('%f, %f' % (expected_envelopes[i][j] - envelope[j], precision / 2))
                return 'fail'
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

    return 'success'

def test_gdal_rasterize_lib_3():

    import test_cli_utilities
    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' ../gdrivers/data/n43.dt0 tmp/n43dt0.shp -i 10 -3d')

    with gdaltest.error_handler():
        ds = gdal.Rasterize('/vsimem/bogus.tif', 'tmp/n43dt0.shp')
    assert ds is None, 'did not expected success'

    ds = gdal.Rasterize('', 'tmp/n43dt0.shp', format='MEM', outputType=gdal.GDT_Byte, useZ=True, layers=['n43dt0'], width=121, height=121, noData=0)

    ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource('tmp/n43dt0.shp')

    ds_ref = gdal.Open('../gdrivers/data/n43.dt0')

    assert ds.GetRasterBand(1).GetNoDataValue() == 0.0, \
        'did not get expected nodata value'

    assert ds.RasterXSize == 121 and ds.RasterYSize == 121, \
        'did not get expected dimensions'

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    for i in range(6):
        assert abs(gt[i] - gt_ref[i]) <= 1e-6, 'did not get expected geotransform'

    wkt = ds.GetProjectionRef()
    assert wkt.find("WGS_1984") != -1, 'did not get expected SRS'

def test_gdal_rasterize_3():

    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    if test_cli_utilities.get_gdal_rasterize_path() is None:
        pytest.skip()

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' ../gdrivers/data/n43.dt0 tmp/n43dt0.shp -i 10 -3d')

    gdaltest.runexternal(test_cli_utilities.get_gdal_rasterize_path() + ' -3d tmp/n43dt0.shp tmp/n43dt0.tif -l n43dt0 -ts 121 121 -a_nodata 0 -q')

    ds_ref = gdal.Open('../gdrivers/data/n43.dt0')
    ds = gdal.Open('tmp/n43dt0.tif')

    assert ds.GetRasterBand(1).GetNoDataValue() == 0.0, \
        'did not get expected nodata value'

    assert ds.RasterXSize == 121 and ds.RasterYSize == 121, \
        'did not get expected dimensions'

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    for i in range(6):
        assert abs(gt[i] - gt_ref[i]) <= 1e-6, 'did not get expected geotransform'

    wkt = ds.GetProjectionRef()
    assert wkt.find("WGS_1984") != -1, 'did not get expected SRS'

def test_gdal_rasterize_4():

    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    if test_cli_utilities.get_gdal_rasterize_path() is None:
        pytest.skip()

    gdal.GetDriverByName('GTiff').Delete('tmp/n43dt0.tif')

    gdaltest.runexternal(test_cli_utilities.get_gdal_rasterize_path() + ' -3d tmp/n43dt0.shp tmp/n43dt0.tif -l n43dt0 -tr 0.008333333333333  0.008333333333333 -a_nodata 0 -a_srs EPSG:4326')

    ds_ref = gdal.Open('../gdrivers/data/n43.dt0')
    ds = gdal.Open('tmp/n43dt0.tif')

    assert ds.GetRasterBand(1).GetNoDataValue() == 0.0, \
        'did not get expected nodata value'

    # Allow output to grow by 1/2 cell, as per #6058
    assert ds.RasterXSize == 122 and ds.RasterYSize == 122, \
        'did not get expected dimensions'

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    assert abs(gt[1] - gt_ref[1]) <= 1e-6 and abs(gt[5] - gt_ref[5]) <= 1e-6, \
        'did not get expected geotransform(dx/dy)'

    # Allow output to grow by 1/2 cell, as per #6058
    assert (abs(gt[0] + (gt[1] / 2) - gt_ref[0]) <= 1e-6 and \
       abs(gt[3] + (gt[5] / 2) - gt_ref[3]) <= 1e-6), \
        'did not get expected geotransform'

    wkt = ds.GetProjectionRef()
    assert wkt.find("WGS_1984") != -1, 'did not get expected SRS'

def test_gdal_contour_5():
    if test_cli_utilities.get_gdal_contour_path() is None:
        return 'skip'

    ds = None

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 10 data/contour_orientation.tif tmp/contour_orientation2.shp')

    ds = ogr.Open('tmp/contour_orientation2.shp')

    expected_contours = [ 'LINESTRING (0 2,'+
                                      '0.5 2.0,'+
                                      '1.5 2.0,'+
                                      '1.954542932445554 2.5,'+
                                      '2.124997615823304 3.5,'+
                                      '1.5 3.954546085074803,'+
                                      '0.5 4.066665649414062,'+
                                      '0.0 4.066665649414062)' ]
    expected_elev = [ 140 ]

    lyr = ds.ExecuteSQL("select * from contour_orientation2 order by elev asc")

    if lyr.GetFeatureCount() != len(expected_contours):
        print('Got %d features. Expected %d' % (lyr.GetFeatureCount(), len(expected_contours)))
        return 'fail'

    i = 0
    test_failed = False
    feat = lyr.GetNextFeature()
    while feat is not None:
        expected_geom = ogr.CreateGeometryFromWkt(expected_contours[i])
        if feat.GetField('elev') != expected_elev[i]:
            print('Got %f. Expected %f' % (feat.GetField('elev'), expected_elev[i]))
            return 'fail'
        if ogrtest.check_feature_geometry(feat, expected_geom) != 0:
            print('Got      %s.\nExpected %s' % (feat.GetGeometryRef().ExportToWkt(),expected_contours[i]))
            test_failed = True
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

    if test_failed:
        return 'fail'
    else:
        return 'success'

def test_gdal_rasterize_lib_3():

    import test_cli_utilities
    if test_cli_utilities.get_gdal_contour_path() is None:
        return 'skip'

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' ../gdrivers/data/n43.dt0 tmp/n43dt0.shp -i 10 -3d')

    with gdaltest.error_handler():
        ds = gdal.Rasterize('/vsimem/bogus.tif', 'tmp/n43dt0.shp')
    if ds is not None:
        gdaltest.post_reason('did not expected success')
        return 'fail'

    ds = gdal.Rasterize('', 'tmp/n43dt0.shp', format = 'MEM', outputType = gdal.GDT_Byte, useZ = True, layers = ['n43dt0'], width = 121, height = 121, noData = 0)

    ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource( 'tmp/n43dt0.shp' )

    ds_ref = gdal.Open('../gdrivers/data/n43.dt0')

    if ds.GetRasterBand(1).GetNoDataValue() != 0.0:
        gdaltest.post_reason('did not get expected nodata value')
        return 'fail'

    if ds.RasterXSize != 121 or ds.RasterYSize != 121:
        gdaltest.post_reason('did not get expected dimensions')
        return 'fail'

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    for i in range(6):
        if (abs(gt[i]-gt_ref[i])>1e-6):
            gdaltest.post_reason('did not get expected geotransform')
            print(gt)
            print(gt_ref)
            return 'fail'

    wkt = ds.GetProjectionRef()
    if wkt.find("WGS_1984") == -1:
        gdaltest.post_reason('did not get expected SRS')
        print(wkt)
        return 'fail'

    return 'success'

def test_gdal_rasterize_3():

    if test_cli_utilities.get_gdal_contour_path() is None:
        return "skip"

    if test_cli_utilities.get_gdal_rasterize_path() is None:
        return "skip"

    gdaltest.runexternal(
        test_cli_utilities.get_gdal_contour_path() + " ../gdrivers/data/n43.dt0 tmp/n43dt0.shp -i 10 -3d"
    )

    gdaltest.runexternal(
        test_cli_utilities.get_gdal_rasterize_path()
        + " -3d tmp/n43dt0.shp tmp/n43dt0.tif -l n43dt0 -ts 121 121 -a_nodata 0 -q"
    )

    ds_ref = gdal.Open("../gdrivers/data/n43.dt0")
    ds = gdal.Open("tmp/n43dt0.tif")

    if ds.GetRasterBand(1).GetNoDataValue() != 0.0:
        gdaltest.post_reason("did not get expected nodata value")
        return "fail"

    if ds.RasterXSize != 121 or ds.RasterYSize != 121:
        gdaltest.post_reason("did not get expected dimensions")
        return "fail"

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    for i in range(6):
        if abs(gt[i] - gt_ref[i]) > 1e-6:
            gdaltest.post_reason("did not get expected geotransform")
            print(gt)
            print(gt_ref)
            return "fail"

    wkt = ds.GetProjectionRef()
    if wkt.find("WGS_1984") == -1:
        gdaltest.post_reason("did not get expected SRS")
        print(wkt)
        return "fail"

    return "success"

def test_gdal_contour_3():
    if test_cli_utilities.get_gdal_contour_path() is None:
        return 'skip'

    try:
        os.remove('tmp/contour.shp')
    except:
        pass
    try:
        os.remove('tmp/contour.dbf')
    except:
        pass
    try:
        os.remove('tmp/contour.shx')
    except:
        pass

    # put -3d just after -fl to test #2793
    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 50 ../gdrivers/data/n43.dt0 tmp/contour.shp')

    ds = ogr.Open('tmp/contour.shp')

    lyr = ds.ExecuteSQL("select distinct elev from contour order by elev asc")

    expected_heights = [ 100, 150, 200, 250, 300, 350, 400, 450 ]
    if lyr.GetFeatureCount() != len(expected_heights):
        print('Got %d features. Expected %d' % (lyr.GetFeatureCount(), len(expected_heights)))
        return 'fail'

    i = 0
    feat = lyr.GetNextFeature()
    while feat is not None:
        if feat.GetField('elev') != expected_heights[i]:
            return 'fail'
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

    return 'success'

def test_gdal_contour_cleanup():
    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource('tmp/contour.shp')
    ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource('tmp/contour_orientation1.shp')
    ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource('tmp/contour_orientation2.shp')
    try:
        os.remove('tmp/gdal_contour.tif')
        os.remove('tmp/gdal_contour_orientation.tif')
    except OSError:
        pass

def test_gdal_contour_5():
    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    ds = None

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 10 data/contour_orientation.tif tmp/contour_orientation2.shp')

    ds = ogr.Open('tmp/contour_orientation2.shp')

    expected_contours =     ['LINESTRING (0.0 1.999999,' +
                             '0.5 1.999999,' +
                             '1.5 1.999999,' +
                             '1.95454293244555 2.5,' +
                             '2.1249976158233 3.5,' +
                             '1.5 3.9545460850748,' +
                             '0.5 4.06666564941406,' +
                             '0.0 4.06666564941406)']
    expected_elev = [140]

    lyr = ds.ExecuteSQL("select * from contour_orientation2 order by elev asc")

    assert lyr.GetFeatureCount() == len(expected_contours)

    i = 0
    test_failed = False
    feat = lyr.GetNextFeature()
    while feat is not None:
        expected_geom = ogr.CreateGeometryFromWkt(expected_contours[i])
        assert feat.GetField('elev') == expected_elev[i]
        if ogrtest.check_feature_geometry(feat, expected_geom) != 0:
            print('Got      %s.\nExpected %s' % (feat.GetGeometryRef().ExportToWkt(), expected_contours[i]))
            test_failed = True
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

    assert not test_failed

def test_gdal_rasterize_3():

    if test_cli_utilities.get_gdal_contour_path() is None:
        return 'skip'

    if test_cli_utilities.get_gdal_rasterize_path() is None:
        return 'skip'

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' ../gdrivers/data/n43.dt0 tmp/n43dt0.shp -i 10 -3d')

    gdaltest.runexternal(test_cli_utilities.get_gdal_rasterize_path() + ' -3d tmp/n43dt0.shp tmp/n43dt0.tif -l n43dt0 -ts 121 121 -a_nodata 0 -q')

    ds_ref = gdal.Open('../gdrivers/data/n43.dt0')
    ds = gdal.Open('tmp/n43dt0.tif')

    if ds.GetRasterBand(1).GetNoDataValue() != 0.0:
        gdaltest.post_reason('did not get expected nodata value')
        return 'fail'

    if ds.RasterXSize != 121 or ds.RasterYSize != 121:
        gdaltest.post_reason('did not get expected dimensions')
        return 'fail'

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    for i in range(6):
        if (abs(gt[i]-gt_ref[i])>1e-6):
            gdaltest.post_reason('did not get expected geotransform')
            print(gt)
            print(gt_ref)
            return 'fail'

    wkt = ds.GetProjectionRef()
    if wkt.find("WGS_1984") == -1:
        gdaltest.post_reason('did not get expected SRS')
        print(wkt)
        return 'fail'

    return 'success'

def test_gdal_rasterize_4():

    if test_cli_utilities.get_gdal_contour_path() is None:
        return "skip"

    if test_cli_utilities.get_gdal_rasterize_path() is None:
        return "skip"

    gdal.GetDriverByName("GTiff").Delete("tmp/n43dt0.tif")

    gdaltest.runexternal(
        test_cli_utilities.get_gdal_rasterize_path()
        + " -3d tmp/n43dt0.shp tmp/n43dt0.tif -l n43dt0 -tr 0.008333333333333  0.008333333333333 -a_nodata 0 -a_srs EPSG:4326"
    )

    ds_ref = gdal.Open("../gdrivers/data/n43.dt0")
    ds = gdal.Open("tmp/n43dt0.tif")

    if ds.GetRasterBand(1).GetNoDataValue() != 0.0:
        gdaltest.post_reason("did not get expected nodata value")
        return "fail"

    # Allow output to grow by 1/2 cell, as per #6058
    if ds.RasterXSize != 122 or ds.RasterYSize != 122:
        gdaltest.post_reason("did not get expected dimensions")
        return "fail"

    gt_ref = ds_ref.GetGeoTransform()
    gt = ds.GetGeoTransform()
    if abs(gt[1] - gt_ref[1]) > 1e-6 or abs(gt[5] - gt_ref[5]) > 1e-6:
        gdaltest.post_reason("did not get expected geotransform(dx/dy)")
        print(gt)
        print(gt_ref)
        return "fail"

    # Allow output to grow by 1/2 cell, as per #6058
    if abs(gt[0] + (gt[1] / 2) - gt_ref[0]) > 1e-6 or abs(gt[3] + (gt[5] / 2) - gt_ref[3]) > 1e-6:
        gdaltest.post_reason("did not get expected geotransform")
        print(gt)
        print(gt_ref)
        return "fail"

    wkt = ds.GetProjectionRef()
    if wkt.find("WGS_1984") == -1:
        gdaltest.post_reason("did not get expected SRS")
        print(wkt)
        return "fail"

    return "success"

def test_gdal_contour_3():
    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    try:
        os.remove('tmp/contour.shp')
    except OSError:
        pass
    try:
        os.remove('tmp/contour.dbf')
    except OSError:
        pass
    try:
        os.remove('tmp/contour.shx')
    except OSError:
        pass

    # put -3d just after -fl to test #2793
    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 50 ../gdrivers/data/n43.dt0 tmp/contour.shp')

    ds = ogr.Open('tmp/contour.shp')

    lyr = ds.ExecuteSQL("select distinct elev from contour order by elev asc")

    expected_heights = [100, 150, 200, 250, 300, 350, 400, 450]
    assert lyr.GetFeatureCount() == len(expected_heights)

    i = 0
    feat = lyr.GetNextFeature()
    while feat is not None:
        assert feat.GetField('elev') == expected_heights[i]
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

def test_gdal_contour_1():
    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    try:
        os.remove('tmp/contour.shp')
    except OSError:
        pass
    try:
        os.remove('tmp/contour.dbf')
    except OSError:
        pass
    try:
        os.remove('tmp/contour.shx')
    except OSError:
        pass

    drv = gdal.GetDriverByName('GTiff')
    sr = osr.SpatialReference()
    sr.ImportFromEPSG(4326)
    wkt = sr.ExportToWkt()

    size = 160
    precision = 1. / size

    ds = drv.Create('tmp/gdal_contour.tif', size, size, 1)
    ds.SetProjection(wkt)
    ds.SetGeoTransform([1, precision, 0, 50, 0, -precision])

    raw_data = array.array('h', [10 for i in range(int(size / 2))]).tostring()
    for i in range(int(size / 2)):
        ds.WriteRaster(int(size / 4), i + int(size / 4), int(size / 2), 1, raw_data,
                       buf_type=gdal.GDT_Int16,
                       band_list=[1])

    raw_data = array.array('h', [20 for i in range(int(size / 2))]).tostring()
    for i in range(int(size / 4)):
        ds.WriteRaster(int(size / 4) + int(size / 8), i + int(size / 4) + int(size / 8), int(size / 4), 1, raw_data,
                       buf_type=gdal.GDT_Int16,
                       band_list=[1])

    raw_data = array.array('h', [25 for i in range(int(size / 4))]).tostring()
    for i in range(int(size / 8)):
        ds.WriteRaster(int(size / 4) + int(size / 8) + int(size / 16), i + int(size / 4) + int(size / 8) + int(size / 16), int(size / 8), 1, raw_data,
                       buf_type=gdal.GDT_Int16,
                       band_list=[1])

    ds = None

    (_, err) = gdaltest.runexternal_out_and_err(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 10 tmp/gdal_contour.tif tmp/contour.shp')
    assert (err is None or err == ''), 'got error/warning'

    ds = ogr.Open('tmp/contour.shp')

    expected_envelopes = [[1.25, 1.75, 49.25, 49.75],
                          [1.25 + 0.125, 1.75 - 0.125, 49.25 + 0.125, 49.75 - 0.125]]
    expected_height = [10, 20]

    lyr = ds.ExecuteSQL("select * from contour order by elev asc")

    assert lyr.GetSpatialRef().ExportToWkt() == wkt, 'Did not get expected spatial ref'

    assert lyr.GetFeatureCount() == len(expected_envelopes)

    i = 0
    feat = lyr.GetNextFeature()
    while feat is not None:
        envelope = feat.GetGeometryRef().GetEnvelope()
        assert feat.GetField('elev') == expected_height[i]
        for j in range(4):
            if abs(expected_envelopes[i][j] - envelope[j]) > precision / 2 * 1.001:
                print('i=%d, wkt=%s' % (i, feat.GetGeometryRef().ExportToWkt()))
                print(feat.GetGeometryRef().GetEnvelope())
                pytest.fail('%f, %f' % (expected_envelopes[i][j] - envelope[j], precision / 2))
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

def test_gdal_contour_4():
    if test_cli_utilities.get_gdal_contour_path() is None:
        pytest.skip()

    try:
        os.remove('tmp/contour_orientation.shp')
    except OSError:
        pass
    try:
        os.remove('tmp/contour_orientation.dbf')
    except OSError:
        pass
    try:
        os.remove('tmp/contour_orientation.shx')
    except OSError:
        pass

    drv = gdal.GetDriverByName('GTiff')
    wkt = 'GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AUTHORITY[\"EPSG\",\"4326\"]]'

    size = 160
    precision = 1. / size

    ds = drv.Create('tmp/gdal_contour_orientation.tif', size, size, 1)
    ds.SetProjection(wkt)
    ds.SetGeoTransform([1, precision, 0, 50, 0, -precision])

# Make the elevation 15 for the whole image
    raw_data = array.array('h', [15 for i in range(int(size))]).tostring()
    for i in range(int(size)):
        ds.WriteRaster(0, i, int(size), 1, raw_data,
                       buf_type=gdal.GDT_Int16,
                       band_list=[1])

# Create a hill with elevation 25
    raw_data = array.array('h', [25 for i in range(2)]).tostring()
    for i in range(2):
        ds.WriteRaster(int(size / 4) + int(size / 8) - 1, i + int(size / 2) - 1, 2, 1, raw_data,
                       buf_type=gdal.GDT_Int16,
                       band_list=[1])

# Create a depression with elevation 5
    raw_data = array.array('h', [5 for i in range(2)]).tostring()
    for i in range(2):
        ds.WriteRaster(int(size / 2) + int(size / 8) - 1, i + int(size / 2) - 1, 2, 1, raw_data,
                       buf_type=gdal.GDT_Int16,
                       band_list=[1])

    ds = None

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 10 tmp/gdal_contour_orientation.tif tmp/contour_orientation1.shp')

    ds = ogr.Open('tmp/contour_orientation1.shp')

    expected_contours = ['LINESTRING (' +
                         '1.628125 49.493749999999999,' +
                         '1.63125 49.496875000000003,' +
                         '1.63125 49.503124999999997,' +
                         '1.628125 49.50625,' +
                         '1.621875 49.50625,' +
                         '1.61875 49.503124999999997,' +
                         '1.61875 49.496875000000003,' +
                         '1.621875 49.493749999999999,' +
                         '1.628125 49.493749999999999)',
                         'LINESTRING (' +
                         '1.38125 49.496875000000003,' +
                         '1.378125 49.493749999999999,' +
                         '1.371875 49.493749999999999,' +
                         '1.36875 49.496875000000003,' +
                         '1.36875 49.503124999999997,' +
                         '1.371875 49.50625,' +
                         '1.378125 49.50625,' +
                         '1.38125 49.503124999999997,' +
                         '1.38125 49.496875000000003)']
    expected_elev = [10, 20]

    lyr = ds.ExecuteSQL("select * from contour_orientation1 order by elev asc")

    assert lyr.GetFeatureCount() == len(expected_contours)

    i = 0
    test_failed = False
    feat = lyr.GetNextFeature()
    while feat is not None:
        expected_geom = ogr.CreateGeometryFromWkt(expected_contours[i])
        assert feat.GetField('elev') == expected_elev[i]
        if ogrtest.check_feature_geometry(feat, expected_geom, 0.01) != 0:
            print('Got      %s.\nExpected %s' % (feat.GetGeometryRef().ExportToWkt(), expected_contours[i]))
            test_failed = True
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

    assert not test_failed

def test_gdal_contour_1():
    if test_cli_utilities.get_gdal_contour_path() is None:
        return 'skip'

    try:
        os.remove('tmp/contour.shp')
    except:
        pass
    try:
        os.remove('tmp/contour.dbf')
    except:
        pass
    try:
        os.remove('tmp/contour.shx')
    except:
        pass

    drv = gdal.GetDriverByName('GTiff')
    wkt = 'GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AUTHORITY[\"EPSG\",\"4326\"]]'

    size = 160
    precision = 1. / size

    ds = drv.Create('tmp/gdal_contour.tif', size, size, 1)
    ds.SetProjection( wkt )
    ds.SetGeoTransform( [ 1, precision, 0, 50, 0, -precision ] )

    raw_data = array.array('h',[10 for i in range(int(size/2))]).tostring()
    for i in range(int(size/2)):
        ds.WriteRaster( int(size/4), i+int(size/4), int(size/2), 1, raw_data,
                        buf_type = gdal.GDT_Int16,
                        band_list = [1] )

    raw_data = array.array('h',[20 for i in range(int(size/2))]).tostring()
    for i in range(int(size/4)):
        ds.WriteRaster( int(size/4)+int(size/8), i+int(size/4)+int(size/8), int(size/4), 1, raw_data,
                        buf_type = gdal.GDT_Int16,
                        band_list = [1] )

    raw_data = array.array('h',[25 for i in range(int(size/4))]).tostring()
    for i in range(int(size/8)):
        ds.WriteRaster( int(size/4)+int(size/8)+int(size/16), i+int(size/4)+int(size/8)+int(size/16), int(size/8), 1, raw_data,
                        buf_type = gdal.GDT_Int16,
                        band_list = [1] )

    ds = None

    gdaltest.runexternal(test_cli_utilities.get_gdal_contour_path() + ' -a elev -i 10 tmp/gdal_contour.tif tmp/contour.shp')

    ds = ogr.Open('tmp/contour.shp')

    expected_envelopes = [ [ 1.25, 1.75, 49.25, 49.75 ],
                           [ 1.25+0.125, 1.75-0.125, 49.25+0.125, 49.75-0.125 ] ]
    expected_height = [ 10, 20 ]

    lyr = ds.ExecuteSQL("select * from contour order by elev asc")

    if lyr.GetSpatialRef().ExportToWkt().find('GCS_WGS_1984') == -1:
        print('Did not get expected spatial ref')
        return 'fail'

    if lyr.GetFeatureCount() != len(expected_envelopes):
        print('Got %d features. Expected %d' % (lyr.GetFeatureCount(), len(expected_envelopes)))
        return 'fail'

    i = 0
    feat = lyr.GetNextFeature()
    while feat is not None:
        envelope = feat.GetGeometryRef().GetEnvelope()
        if feat.GetField('elev') != expected_height[i]:
            print('Got %f. Expected %f' % (feat.GetField('elev'), expected_height[i]))
            return 'fail'
        for j in range(4):
            if abs(expected_envelopes[i][j] - envelope[j]) > precision/2*1.001:
                print('i=%d, wkt=%s' % (i, feat.GetGeometryRef().ExportToWkt()))
                print(feat.GetGeometryRef().GetEnvelope())
                print(expected_envelopes[i])
                print('%f, %f' % (expected_envelopes[i][j] - envelope[j], precision / 2))
                return 'fail'
        i = i + 1
        feat = lyr.GetNextFeature()

    ds.ReleaseResultSet(lyr)
    ds.Destroy()

    return 'success'