[gdal-dev] comparing two rasters
jdmorgan
jdmorgan at unca.edu
Thu Feb 23 11:32:38 EST 2012
Hi Chaitanya,
This works great. Thanks so much. GDAL/Python is a really powerful
combination.
Cheers,
Derek
On 2/23/2012 9:31 AM, Chaitanya kumar CH wrote:
> 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
> <mailto:jdmorgan at unca.edu>> wrote:
>
> ThanksChaitanya,
>
> 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
>> <mailto: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
>>> <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 <tel:%2B91-9494447584>
>>> 17.2416N 80.1426E
>>
>>
>> --
>> Derek @ NEMAC
>>
>>
>>
>>
>> --
>> Best regards,
>> Chaitanya kumar CH.
>>
>> +91-9494447584 <tel:%2B91-9494447584>
>> 17.2416N 80.1426E
>
>
> --
> Derek @ NEMAC
>
>
>
>
> --
> 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/20120223/4807b748/attachment-0001.html
More information about the gdal-dev
mailing list