Hello GDAL guru's,
I am working on a python script where I read in two rasters of similar
extent and resolution.Then I re-assign any values that are greater that
zero to a 1.Next, I compare to the rasters and attempt to create a third
resulting raster which has 1's everywhere that the two input rasters
match up.However, my results are not exactly as expected. The third
raster only has the values of the second raster. Any help would be
greatly appreciated.Here is the script as it is now:
#!/usr/bin/python
from osgeo import gdal, osr, gdalconst
import os, sys, time
import struct
import numpy as np
np.set_printoptions(threshold='nan')
gdal.AllRegister()
print 'Raster Source 1------------------'
ds1 = gdal.Open('C:\Data\TE300by300.img', gdal.GA_ReadOnly)
cols1 = ds1.RasterXSize
rows1 = ds1.RasterYSize
bands1 = ds1.RasterCount
print "Columns: " + str(cols1)
print "Rows: " + str(rows1)
print "Bands: " + str(bands1)
gt1 = ds1.GetGeoTransform()
width = ds1.RasterXSize
height = ds1.RasterYSize
minx = gt1[0]
print "Left(minx): "+ str(minx)
miny = gt1[3] + width*gt1[4] + height*gt1[5]
print "Bottom(miny): "+ str(miny)
maxx = gt1[0] + width*gt1[1] + height*gt1[2]
print "Right(maxx): "+ str(maxx)
maxy = gt1[3]
print "Top(maxy): "+ str(maxy)
pixWidth = gt1[1]
print "Pixel Width " + str(pixWidth)
pixHeight = gt1[5]
print "Pixel Height " + str(pixHeight)
xOrigin = gt1[0]
yOrigin = gt1[3]
print 'Raster Source 2------------------'
ds2 = gdal.Open('C:\Data\LowElev300by300.img', gdal.GA_ReadOnly)
cols2 = ds2.RasterXSize
rows2 = ds2.RasterYSize
bands2 = ds2.RasterCount
print "Columns: " + str(cols2)
print "Rows: " + str(rows2)
print "Bands: " + str(bands2)
gt2 = ds2.GetGeoTransform()
width = ds2.RasterXSize
height = ds2.RasterYSize
minx = gt2[0]
print "Left(minx): "+ str(minx)
miny = gt2[3] + width*gt2[4] + height*gt2[5]
print "Bottom(miny): "+ str(miny)
maxx = gt2[0] + width*gt2[1] + height*gt2[2]
print "Right(maxx): "+ str(maxx)
maxy = gt2[3]
print "Top(maxy): "+ str(maxy)
pixWidth = gt2[1]
print "Pixel Width " + str(pixWidth)
pixHeight = gt2[5]
print "Pixel Height " + str(pixHeight)
xOrigin = gt2[0]
yOrigin = gt2[3]
format = "HFA"
dst_file = "out.img"
dst_driver = gdal.GetDriverByName(format)
dst_ds = dst_driver.Create(dst_file, width, height, 1, gdal.GDT_Byte )
#driver.Create( outfile, outwidth, outheight, numbands, gdaldatatype)
empty = np.empty([height, width], np.uint8 )
srs = osr.SpatialReference()
dst_ds.SetProjection(ds2.GetProjection())
dst_ds.SetGeoTransform(ds2.GetGeoTransform())
dst_ds.GetRasterBand(1).WriteArray(empty)
#This is a way to go through a given raster band
#one-by-one as an array by row, getting all of the columns for
for r in range(rows1):
data1 = ds1.GetRasterBand(1).ReadAsArray(0, r, cols1, 1)
data1[data1>0]=1
print "data1: " + str(data1)
data2 = ds2.GetRasterBand(1).ReadAsArray(0, r, cols2, 1)
data2[data2>0]=1
print "data2: " + str(data2)
result_bools = np.logical_and(np.logical_and(data1 != 0, data2 != 0),
data1 == data2)
result_ints = np.array(result_bools, dtype=int)
print "result_ints: " + str(result_ints)
dst_ds.GetRasterBand(1).WriteArray(result_ints, 0, r)
dst_ds = None
Cheers,
Derek
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
