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@xxxxxxxxxxxxxxx
http://lists.osgeo.org/mailman/listinfo/gdal-dev