If you use the terminal (OSX):
Try fixing the GDAL Path.
If you use the GDAL Framework of KyngChaos, the complete Path of the command is /Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize or .../Versions/1.9/...or...
So:
1) in the Terminal:
export PATH=/Library/Frameworks/GDAL.framework/Versions/1.10/Programs
and you can use gdal_rasterizein the terminal during the session.
2) put in the /Users/you/.bash_profile file the line:
export PATH=/Library/Frameworks/GDAL.framework/Versions/1.10/Programs:$PATH
and you can use gdal_rasterize in the terminal all the time
3) in all other cases (script, through IDLE, QGIS, etc.) , you must use the complete Path:
/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize
It is therefore normal that os.system('gdal_rasterize... does not work, you need
os.system('/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize...
If you use Python:
1) os.system is deprecated, you must use the subprocess module which you can control the stdout and stderr in a flexible fashion (look at 17.1.4.3. Replacing os.system(), subprocess.check_output or Calling an external command in Python). You can use subprocess.call(), subprocess.check_call()or subprocess.Popen():
args = ['/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize', '-tr', '7999.8', '9999.53', '-te', '170000.2', '168000.27', '177999.8', '177999.8' '-a Pixelkost /Users/PeterVanvoorden/Desktop/GRBgis_322/Shapefile/Adp322.shp', '/Users/PeterVanvoorden/Desktop/Test.tif']
proc = subprocess.Popen(args, stdout=subprocess.PIPE)
proc.wait() # or proc.communicate()[0] or proc.stdout.read() if you want to read the errors
2) if you use the Python osgeo.gdal module, it has the gdal.RasterizeLayer() function ( Python GDAL/OGR Cookbook,Python: Rasterizing a GDAL layer, Making Rasters From Vectors or GDAL RasterizeLayer Doesn't Burn All Polygons to Raster) and you don't need to call gdal_rasterize.
from osgeo import gdal
raster = gdal.RasterizeLayer(...)
I don't understand your script. Why use PyShp (shapefile), if you use osgeo.ogr ? You could do the same thing with ogr only ( better way to duplicate a layer using ogr in python? ).
If you want to import ogr and gdal from osgeo, the formulation
import osgeo.ogr, osgeo.ogr
from osgeo import ogr
from osgeo import gdal
is a redundancy (you import ogr and gdal twice)
Same for sf.records()and sf.shapeRecords():
sf = shapefile.Reader("strati")
records = sf.records() # Reading the records of the features
for rec in records:
print rec
[30, 130, 'incl']
[55, 145, 'incl']
[40, 155, 'incl']
and
shapRecs = sf.shapeRecords()
for rec in shapRecs:
print rec.record
[30, 130, 'incl']
[55, 145, 'incl']
[40, 155, 'incl']
If you only want the geometry:
for rec in sf.iterShapes():
x,y = rec.points[0]
print x,y
272070.600041, 155389.3879200000
271066.03214800003, 154475.63137700001
273481.498868, 153923.49298800001
and
for rec in shapRecs:
x,y = rec.shape.points[0]
If you just want to copy the shapefile:
w = shapefile.Writer(shapefile.POINT)
sf = shapefile.Reader("strati")
w.fields = list(sf.fields)
for rec in sf.records():
w.records.append(rec)
w._shapes.extend(sf.shapes())
w.save("duplicate")
If you want to add a Field, look at Add a Field to an Existing Shapefile
If you want lists, the easiest way is:
sf = shapefile.Reader("yourfile.shp")
shapes = sf.shapes() # -> list of all the geometries
fields = sf.fields[1:] #-> fields definition in a list (without the Deletion Flag)
field_names = = [field[0] for field in fields] #-> list of fields names
attributes = sf.records() #-> list with all the attributes
Best Answer
You can calculate a buffer directly with
ogr2ogrusing GDAL >= 1.10 with SQLite and SpatiaLite support. If the raster spatial resolution is e.g. 1m:Finally, you can adapt and apply the first example of
gdal_rasterizedocumentation: