1

I am using following code to take subset of all files in a folder. The code is not giving any error but it is not creating output files. How to solve it?

import numpy as np
import re
import gdal
import ogr
from osgeo import gdal_array
from gdalconst import *
import os
import osr
import subprocess
e=1
indir = '/home/vigna/Documents/polynyas/SIC/2017/Sept'
directory=os.fsencode(indir)
for file in os.listdir(directory):    
    file1=os.fsdecode(file)
    if file1.endswith(".tif"):
        ds = gdal.Open(os.path.join(indir,file1),GA_ReadOnly)
        out=str(e)+'_subset.tif'
        inputf = file1# input raster
        #print(len(inDS))
        inDS=str(inputf[0:33])
        #print(inDS)
        outDS = out# output raster
        #print(outDS)
        ulx = -161255
        uly = 2971737
        lrx = 437613
        lry = 2475328
        translate = 'gdal_translate -projwin %s %s %s %s %s %s' %(ulx, uly, lrx, lry, inDS, outDS)
        #print(translate)
        os.system(translate)
        e+=1
  • not the exact solution, but you should try with subprocess.call instead of os.system. And check the message of the "printed" command line + use full path for the output file name. – radouxju Mar 27 '18 at 6:19
  • Try using subprocess as @radouxju indicated. If you do it like suggested here, you should see something, especially when you drop the if and print stdout and sdterr no matter what p.returncode is – s6hebern Mar 29 '18 at 7:54
1

The problem is that you are calling the system's gdal_translate binary from within the python environment and passing to it a python (type = gdal.Dataset) object, which it does not recognize as a valid format (gdalinfo --formats for the supported list). For what you are doing, you do not need to import/use the python bindings, but simply build the command in python and pass it off to the OS. This version of your code should work:

import os,sys
e=1
indir = '/home/vigna/Documents/polynyas/SIC/2017/Sept'
directory=os.fsencode(indir)
for file in os.listdir(directory):
    file1=os.fsdecode(file)
    if file1.endswith(".tif"): # I would check the file type this way: "if os.path.splitext(file1)[1] =='.tif':"
        inputf = file1# input raster
        out=str(e)+'_subset.tif' # I would retain the original filename for clarity i.e. "out = os.path.basename(file1).splitext[0]+'_subset.tif'"
        ulx = -161255
        uly = 2971737
        lrx = 437613
        lry = 2475328
        translate = 'gdal_translate -projwin %s %s %s %s %s %s' %(ulx, uly, lrx, lry, file1, out)
        os.system(translate)
        e+=1
1

I think @planetsandman has you covered, but I wanted to add that you can debug these system calls yourself if you use

import subprocess
translate = 'your cmd line call string'
proc = subprocess.Popen(translate, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
stdout,stderr=proc.communicate()

The stdout variable will contain the output of the call, and the stderr variable will contain any error messages.

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.