Calling Gdal contour from python/ipython

Question

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.

• artwork21: Your dst variable does not look correct, you are appending file path with an extension with a new string with an extension (e.g. C:/Temp/h_Max.ASCb.shp)
• SimonDeS: yes this is poor programming but it does actually work - just produces a file with .ASC in the file name. I have manually ran the command in command line with both the .ASC in and it taken out, produces the same shapefile.
• gisnside: Does this command work on a single asc when you replace by the right values ? > gdalcontour -fl -888 -inodata -f "ESRI Shapefile "' + src + " " + dst [EDIT : you replied just before i send my comment :) ]
• gisnside: Does this work for example (you can replace -i 25 by your favorite value of interval -i 10 or -i 5, etc) : gdal_contour -i 25 -inodata -f "ESRI shapefile" "C:\D\OUT\Wray_0030F_012_h_Max.asc" "C:\D\OUT\Wray_0030F_012_h_Max_25.shp" (if yes, then put it in your code and retry)
• SimonDeS: Microsoft Windows [Version 6.3.9600] (c) 2013 Microsoft Corporation. All rights reserved. C:\Users\Simon.DeSmit>gdal_contour.exe -fl -888 -inodata -f "ESRI Shapefile" C:\ D\OUT\Wray_0010F_012_h_Max.asc C:\D\OUT\Wray_0010F_012_h_Max.asc.shp 0...10...20...30...40...50...60...70...80...90...100 - done. yes. I am using fl -888 rather than -i 25 as i want to draw a line around the entire dataset sans no data.

Answer 1

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.

• anakaine: For posterity: This answer doesn't ouput a projection file, so the resulting shp is not georeferenced.