[gdal-dev] comparing two rasters

[Chaitanya kumar CH]

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 at unca.edu> 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 at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>



-- 
Best regards,
Chaitanya kumar CH.

+91-9494447584
17.2416N 80.1426E
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20120222/0195c5cb/attachment.html