I am trying to use a subprocess.call method in python to use 'gdalwarp'
on a .img file. So far everything that I've put in my script runs
fine...but the gdalwarp. There is no error...so I assume I am somehow
not understanding the gdalwarp statement as well as I thought I might.
Anybody see what I am missing (or can direct me to a strong reference
regarding gdalwarp)?
Enclosing full code in .txt.
subprocess.call('gdalwarp %s %s -t_srs "+proj=lcc
+lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5
+x_0=200000 +y_0=750000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"'
%(raster, 'new_'+raster))
# osgeo package contains the GDAl, OGR, and OSR libraries
from osgeo import gdal
import osr
import os
import fnmatch
import sys
import string
import re
import numpy
import subprocess
# To begin, I set the initial directory workspace. "directory" is the initial
directory to
# search under,
# WARNING : THIS MUST BE SET BEFORE THE USER BEGINS USING THIS SCRIPT
folder = r'C:\WorkOrder\WO_61_HUD_FFEA\hilshd_SKC_10m\img\n42w074'
# Using the os.walk and fnmatch functions, I construct a list of raster files
# to pass in my gdal.open() funtion. At the moment, I've set the file extension
to '.flt'.
# WARNING : THE 'PATTERN' VARIABLE MUST BE SET BEFORE THE USER BEGINS USING
# THIS SCRIPT
pattern = "*.img"
rastList = []
# create a list of files that are within the 'folder' directory
for path, dirs, files in os.walk(folder):
for filename in fnmatch.filter(files, pattern):
rastList.append(os.path.join(path, filename))
# It appears that unlike OGR (which needs the ogr.GetDriverByName())
# in python GDAL detects what type of file you are opening. It loads and
# registers the correct driver, so I can just proceed to using my for loop
for raster in rastList:
datafile = gdal.Open(raster)
# command to extract information on the img size
cols = datafile.RasterXSize
rows = datafile.RasterYSize
bands = datafile.RasterCount
# access the projection information using GetProjection()
# returns info in a WKT string format
oldprojInfo = datafile.GetProjection()
# create spatial refterence object for old projection
spatialRef = osr.SpatialReference()
spatialRef.ImportFromWkt(oldprojInfo)
spatialRefProj = spatialRef.ExportToProj4()
# transform the the old projection using the subprocess.call
# and gdalwarp, could add multiprocessing...see
#
http://gis.stackexchange.com/questions/6669/converting-projected-geotiff-to-wgs84-with-gdal-and-python
subprocess.call('gdalwarp %s %s -t_srs "+proj=lcc
+lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5
+x_0=200000 +y_0=750000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"' %(raster,
'new_'+raster))
# access the projection information using GetProjection()
newprojInfo = datafile.GetProjection()
# create a spatial refrence object for new projection
spatialRefNew = osr.SpatialReference()
spatialRefNew.ImportFromWkt(newprojInfo)
spatialRefNewProj = spatialRefNew.ExportToProj4()
# print information to shell screen
print "Number of columns: " + str(cols)
print "Number of rows: " + str(rows)
print "Number of bands: " + str(bands)
print "The original projection is (in WKT): " + oldprojInfo
print "The original projection is (in Proj4): " + str(spatialRefProj)
print "the new projection is (in WKT): " + newprojInfo
print "the new projection is (in Proj4): " + str(spatialRefNewProj)
if datafile is None:
print 'Could not open ' + raster
sys.exit(1)
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
