[gdal-dev] comparing two rasters

[Chaitanya kumar CH]

Try this Derek,

for r in range(rows1):
    data1 = ds1.GetRasterBand(1).ReadAsArray(0, r, cols1, 1)
    print "data1: " + str(data1)
    data2 = ds2.GetRasterBand(1).ReadAsArray(0, r, cols2, 1)
    print "data2: " + str(data2)
    result_bools = np.logical_and((data1 > 0), (data2 > 0))
    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

On Wed, Feb 22, 2012 at 5:26 PM, jdmorgan <jdmorgan at unca.edu> wrote:

>  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 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
>
>
>
> --
> Derek @ NEMAC
>
>


-- 
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/6ab3a11f/attachment-0001.html