[gdal-dev] comparing two rasters

[jdmorgan]

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 
> <mailto: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 <mailto: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

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