2

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.

  • 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) – artwork21 Sep 13 '16 at 14:12
  • 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. – SimonDeS Sep 13 '16 at 14:14
  • 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 Sep 13 '16 at 14:14
  • 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) – gisnside Sep 13 '16 at 14:22
  • 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. – SimonDeS Sep 13 '16 at 14:50
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.