[gdal-dev] comparing two rasters

Chaitanya kumar CH chaitanya.ch at gmail.com
Thu Feb 23 09:31:29 EST 2012 

Derek,

You can assign a color table.
Ex: (0,0,0) for 0, (255,255,255) for 1.
You can assign a color table to a raster band using the method
SetRasterColorTable().
You can find some sample python code in GDAL's test suite[1][2].

[1]: http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/colortable.py
[2]:
http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/tiff_write.py#L916

On Thu, Feb 23, 2012 at 7:45 PM, jdmorgan <jdmorgan at unca.edu> wrote:

>  Thanks Chaitanya,
>
> Actually this seems to be working, which is pretty cool. I think where I
> am getting confused is that I am attempting to verify the results by
> loading the resulting raster data into a GIS e.g. ArcGIS or QGIS and doing
> spot checking.  However, as I haven’t done any actual color coding on the
> classes e.g. making the 1 black and 0’s white, I think I am not able to
> verify it correctly.  If I let ArcGIS do a default classification it does
> something line 0’s white and 0-1 black.  Is there a way to color class
> this in GDAL python?  That would probably be a better approach.
>
> Cheers,
>
> Derek
>
>
> On 2/22/2012 10:01 AM, Chaitanya kumar CH wrote:
>
> 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
>
>
>
> --
> 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/20120223/e411be00/attachment.html

More information about the gdal-dev mailing list