[gdal-dev] comparing two rasters

[jdmorgan]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20120221/aecb1250/attachment-0001.html