Hi Chaitanya,
I am using data1[data1>0]=1 to convert any of the values in the row of data that is greater than 0 to a one. I am doing this because the values are varied, but I am only interested in the fact that there is a value at all. My end goal is to compare the two input rasters for places where they have any data (not zero). But, as I mentioned I am certainly stuck somewhere getting results that I don't expect.

Thanks,
Derek

On 2/21/2012 11:33 PM, Chaitanya kumar CH wrote:
Derek,

Can you explain the following lines towards the bottom of the script?
data1[data1>0]=1
...
data2[data2>0]=1


On Wed, Feb 22, 2012 at 8:46 AM, jdmorgan <jdmorgan@xxxxxxxx <mailto:jdmorgan@xxxxxxxx>> wrote:

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




--
Best regards,
Chaitanya kumar CH.

+91-9494447584
17.2416N 80.1426E


--
Derek @ NEMAC

_______________________________________________
gdal-dev mailing list
gdal-dev@xxxxxxxxxxxxxxx
http://lists.osgeo.org/mailman/listinfo/gdal-dev