[gdal-dev] comparing two rasters
jdmorgan
jdmorgan at unca.edu
Thu Feb 23 09:15:48 EST 2012
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
> 17.2416N 80.1426E
--
Derek @ NEMAC
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20120223/32233c4c/attachment-0001.html
More information about the gdal-dev
mailing list