Guest User

Untitled

a guest
Oct 26th, 2016
54
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.11 KB | None | 0 0
  1. """
  2. Example of raster reclassification using OpenSource Geo Python
  3.  
  4. """
  5. import numpy, sys
  6. from osgeo import gdal
  7. from osgeo.gdalconst import *
  8.  
  9.  
  10. # register all of the GDAL drivers
  11. gdal.AllRegister()
  12.  
  13. # open the image
  14. inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
  15. if inDs is None:
  16. print 'Could not open image file'
  17. sys.exit(1)
  18.  
  19. # read in the crop data and get info about it
  20. band1 = inDs.GetRasterBand(1)
  21. rows = inDs.RasterYSize
  22. cols = inDs.RasterXSize
  23.  
  24. cropData = band1.ReadAsArray(0,0,cols,rows)
  25.  
  26. listAg = [1,5,6,22,23,24,41,42,28,37]
  27. listNotAg = [111,195,141,181,121,122,190,62]
  28.  
  29. # create the output image
  30. driver = inDs.GetDriver()
  31. #print driver
  32. outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
  33. if outDs is None:
  34. print 'Could not create reclass_40.tif'
  35. sys.exit(1)
  36.  
  37. outBand = outDs.GetRasterBand(1)
  38. outData = numpy.zeros((rows,cols), numpy.int16)
  39.  
  40.  
  41. for i in range(0, rows):
  42. for j in range(0, cols):
  43.  
  44. if cropData[i,j] in listAg:
  45. outData[i,j] = 100
  46. elif cropData[i,j] in listNotAg:
  47. outData[i,j] = -100
  48. else:
  49. outData[i,j] = 0
  50.  
  51.  
  52. # write the data
  53. outBand.WriteArray(outData, 0, 0)
  54.  
  55. # flush data to disk, set the NoData value and calculate stats
  56. outBand.FlushCache()
  57. outBand.SetNoDataValue(-99)
  58.  
  59. # georeference the image and set the projection
  60. outDs.SetGeoTransform(inDs.GetGeoTransform())
  61. outDs.SetProjection(inDs.GetProjection())
  62.  
  63. del outData
  64.  
  65. #import arcgisscripting
  66.  
  67. gp = arcgisscripting.create(9.3)
  68.  
  69. gp.AddToolbox("SA") # addint spatial analyst toolbox
  70.  
  71. rasterA = @"C:rasterA.tif"
  72. rasterB = @"C:rasterB.tif"
  73.  
  74. rasterC = @"C:rasterC.tif" # this raster does not yet exist
  75. rasterD = @"C:rasterD.tif" # this raster does not yet exist
  76.  
  77. gp.Minus_SA(rasterA,rasterB,rasterC)
  78.  
  79. gp.Times_SA(rasterA,rasterB,rasterD)
  80.  
  81. # lets try to use more complex functions
  82.  
  83. # lets build and expression first
  84.  
  85. expression1 = "slope( " + rasterC + ")"
  86. expression2 = "(" + rasterC " + " rasterD + ") - " + rasterA
  87.  
  88. gp.SingleOutputMapAlgebra_SA(expression1,@"C:result_exp1.tif")
  89. gp.SingleOutputMapAlgebra_SA(expression2,@"C:result_exp2.tif")
Add Comment
Please, Sign In to add comment