Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- """
- Example of raster reclassification using OpenSource Geo Python
- """
- import numpy, sys
- from osgeo import gdal
- from osgeo.gdalconst import *
- # register all of the GDAL drivers
- gdal.AllRegister()
- # open the image
- inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
- if inDs is None:
- print 'Could not open image file'
- sys.exit(1)
- # read in the crop data and get info about it
- band1 = inDs.GetRasterBand(1)
- rows = inDs.RasterYSize
- cols = inDs.RasterXSize
- cropData = band1.ReadAsArray(0,0,cols,rows)
- listAg = [1,5,6,22,23,24,41,42,28,37]
- listNotAg = [111,195,141,181,121,122,190,62]
- # create the output image
- driver = inDs.GetDriver()
- #print driver
- outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
- if outDs is None:
- print 'Could not create reclass_40.tif'
- sys.exit(1)
- outBand = outDs.GetRasterBand(1)
- outData = numpy.zeros((rows,cols), numpy.int16)
- for i in range(0, rows):
- for j in range(0, cols):
- if cropData[i,j] in listAg:
- outData[i,j] = 100
- elif cropData[i,j] in listNotAg:
- outData[i,j] = -100
- else:
- outData[i,j] = 0
- # write the data
- outBand.WriteArray(outData, 0, 0)
- # flush data to disk, set the NoData value and calculate stats
- outBand.FlushCache()
- outBand.SetNoDataValue(-99)
- # georeference the image and set the projection
- outDs.SetGeoTransform(inDs.GetGeoTransform())
- outDs.SetProjection(inDs.GetProjection())
- del outData
- #import arcgisscripting
- gp = arcgisscripting.create(9.3)
- gp.AddToolbox("SA") # addint spatial analyst toolbox
- rasterA = @"C:rasterA.tif"
- rasterB = @"C:rasterB.tif"
- rasterC = @"C:rasterC.tif" # this raster does not yet exist
- rasterD = @"C:rasterD.tif" # this raster does not yet exist
- gp.Minus_SA(rasterA,rasterB,rasterC)
- gp.Times_SA(rasterA,rasterB,rasterD)
- # lets try to use more complex functions
- # lets build and expression first
- expression1 = "slope( " + rasterC + ")"
- expression2 = "(" + rasterC " + " rasterD + ") - " + rasterA
- gp.SingleOutputMapAlgebra_SA(expression1,@"C:result_exp1.tif")
- gp.SingleOutputMapAlgebra_SA(expression2,@"C:result_exp2.tif")
Add Comment
Please, Sign In to add comment