[GIS] Calling Gdal contour from python/ipython

gdalgdal-contourpythonpython-2.7

Hello I am trying to patch process a folder full of ASCs to contours using Gdal from python. here is my code;

import os 
import glob

PATH1 = raw_input("folder containing ASCs path please:") 


inputlistB = []

    
for filename in glob.glob(os.path.join(PATH1, '*h_Max.ASC')):
    print filename
    inputlistB.append(filename)
    
print inputlistB

for filepath in inputlistB:
    src = filepath
    dst = filepath+'b.shp'
    os.system('gdalcontour -fl -888 -inodata -f "ESRI Shapefile "' + src + " " + dst)

The ipython console returns this;

folder containing ASCs path please:C:\D\OUT

C:\D\OUT\Wray_0010F_012_h_Max.asc C:\D\OUT\Wray_0020F_012_h_Max.asc

C:\D\OUT\Wray_0030F_012_h_Max.asc

['C:\D\OUT\Wray_0010F_012_h_Max.asc',

'C:\D\OUT\Wray_0020F_012_h_Max.asc',

'C:\D\OUT\Wray_0030F_012_h_Max.asc']

showing that the inputlist functions as intended. However no new gis layers are generated in my target folder. If I take the statement ran by os.system and copy it with file paths into command line, it works. Though it does crash after producing the shapefile. If I just run the final line of the for loop os.system(… the Ipython console simply returns :1

Any idea why gdal would crash after producing a shapefile / successfully running?

Any ideas as to why my code does not produce any results despite no python errors?

Any alternative / better solutions to this problem?

also I use the file path input and loop in a lot of programs, feel free to steal it if this is news to you.

Best Answer

import os
import glob
from osgeo import gdal
from osgeo.gdalconst import *
from numpy import *
from osgeo import ogr


PATH1 = raw_input("folder containing ASCs path please:") 


inputlistB = []


for filename in glob.glob(os.path.join(PATH1, '*h_Max.ASC')):
    print filename
    inputlistB.append(filename)



# main loop - open each asci and spit out a contour file at -888level

for ascipath in inputlistB:
    asci = str(ascipath)
    indataset1 = gdal.Open( asci, GA_ReadOnly)
    in1 = indataset1.GetRasterBand(1)
    dst_filename = asci+'contour.shp'
    #Generate layer to save Contourlines in
    ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(dst_filename)
    contour_shp = ogr_ds.CreateLayer('contour')

    field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)
    contour_shp.CreateField(field_defn)
    field_defn = ogr.FieldDefn("elev", ogr.OFTReal)
    contour_shp.CreateField(field_defn)

    #Generate Contourlines
    gdal.ContourGenerate(in1, 0, 0, [-888], 0, 0, contour_shp, 0, 1)
    ogr_ds.Destroy()

Solved. The above code works. It would appear the problem with my original code was in part a pathing issue but I have gone with the suggestions in the comments and cribbed some code from other threads on stackexchange to use gdal within python. Regards to user Ben Mewes for the gdal code. I have modified it so it produces a single outline around my elevation data. It needs a little tlc but this will work for patch processing NCFDD outlines.

Related Question