<html>
<head>
<meta content="text/html; charset=UTF-8" http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Hi Chaitanya, <br>
This works great. Thanks so much. GDAL/Python is a really powerful
combination.<br>
<br>
Cheers, <br>
Derek<br>
<br>
On 2/23/2012 9:31 AM, Chaitanya kumar CH wrote:
<blockquote
cite="mid:CAMKgpOYcBURVemCiNrPjSBRc6_XBWTRSwYinr06P9x0hK-Pnfw@mail.gmail.com"
type="cite">Derek,<br>
<br>
You can assign a color table.<br>
Ex: (0,0,0) for 0, (255,255,255) for 1.<br>
You can assign a color table to a raster band using the method
SetRasterColorTable().<br>
You can find some sample python code in GDAL's test suite[1][2].<br>
<br>
[1]: <a moz-do-not-send="true"
href="http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/colortable.py">http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/colortable.py</a><br>
[2]: <a moz-do-not-send="true"
href="http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/tiff_write.py#L916">http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/tiff_write.py#L916</a><br>
<br>
<div class="gmail_quote">On Thu, Feb 23, 2012 at 7:45 PM, jdmorgan
<span dir="ltr"><<a moz-do-not-send="true"
href="mailto:jdmorgan@unca.edu">jdmorgan@unca.edu</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0
.8ex;border-left:1px #ccc solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000">
<p class="MsoNormal">Thanks<big> </big>Chaitanya, <br>
</p>
<p class="MsoNormal">Actually this seems to be working,
which is pretty cool.<span> </span>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.<span> </span>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.<span> </span>If I let
ArcGIS do a default classification it does something line
0’s white and 0-1 black.<span> </span>Is there a way to
color class this in GDAL python?<span> </span>That would
probably be a better approach. </p>
<p class="MsoNormal">Cheers, </p>
<p class="MsoNormal">Derek</p>
<div>
<div class="h5"> <br>
<br>
On 2/22/2012 10:01 AM, Chaitanya kumar CH wrote:
<blockquote type="cite">Try this Derek,<br>
<br>
<span style="font-family:courier new,monospace">for r
in range(rows1):</span><br
style="font-family:courier new,monospace">
<span style="font-family:courier new,monospace">
data1 = ds1.GetRasterBand(1).ReadAsArray(0, r,
cols1, 1)</span><br style="font-family:courier
new,monospace">
<span style="font-family:courier new,monospace">
print "data1: " + str(data1)</span><br
style="font-family:courier new,monospace">
<span style="font-family:courier new,monospace">
data2 = ds2.GetRasterBand(1).ReadAsArray(0, r,
cols2, 1)</span><br style="font-family:courier
new,monospace">
<span style="font-family:courier new,monospace">
print "data2: " + str(data2)</span><br
style="font-family:courier new,monospace">
<span style="font-family:courier new,monospace">
result_bools = np.logical_and((data1 > 0), (data2
> 0)</span>)<br style="font-family:courier
new,monospace">
<span style="font-family:courier new,monospace">
result_ints = np.array(result_bools, dtype=int)</span><br
style="font-family:courier new,monospace">
<span style="font-family:courier new,monospace">
print "result_ints: " + str(result_ints)</span><br
style="font-family:courier new,monospace">
<span style="font-family:courier new,monospace">
dst_ds.GetRasterBand(1).WriteArray(result_ints, 0,
r)</span><br style="font-family:courier
new,monospace">
<span style="font-family:courier new,monospace">dst_ds
= None</span><br>
<br>
<div class="gmail_quote">On Wed, Feb 22, 2012 at 5:26
PM, jdmorgan <span dir="ltr"><<a
moz-do-not-send="true"
href="mailto:jdmorgan@unca.edu" target="_blank">jdmorgan@unca.edu</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0
.8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF"> Hi
Chaitanya, <br>
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.<br>
<br>
Thanks, <br>
Derek
<div>
<div><br>
<br>
On 2/21/2012 11:33 PM, Chaitanya kumar CH
wrote:
<blockquote type="cite">Derek,<br>
<br>
Can you explain the following lines
towards the bottom of the script?<br>
<span style="font-family:courier
new,monospace">data1</span><span
style="font-family:courier
new,monospace">[</span><span
style="font-family:courier
new,monospace">data1></span><span
style="font-family:courier
new,monospace">0</span><span
style="font-family:courier
new,monospace">]</span><span
style="font-family:courier
new,monospace">=</span><span
style="font-family:courier
new,monospace">1<br>
...<br>
</span><span style="font-family:courier
new,monospace">data2</span><span
style="font-family:courier
new,monospace">[</span><span
style="font-family:courier
new,monospace">data2></span><span
style="font-family:courier
new,monospace">0</span><span
style="font-family:courier
new,monospace">]</span><span
style="font-family:courier
new,monospace">=</span><span
style="font-family:courier
new,monospace">1</span><br>
<br>
<br>
<div class="gmail_quote">On Wed, Feb 22,
2012 at 8:46 AM, jdmorgan <span
dir="ltr"><<a
moz-do-not-send="true"
href="mailto:jdmorgan@unca.edu"
target="_blank">jdmorgan@unca.edu</a>></span>
wrote:<br>
<blockquote class="gmail_quote"
style="margin:0 0 0
.8ex;border-left:1px #ccc
solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000">
<p class="MsoNormal">Hello GDAL
guru’s, </p>
<p class="MsoNormal">I am working on
a python script where I read in
two rasters of similar extent and
resolution.<span> </span>Then I
re-assign any values that are
greater that zero to a 1.<span> </span>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.<span> </span>However,
my results are not exactly as
expected.<span> The third raster
only has the values of the
second raster. </span>Any help
would be greatly appreciated.<span>
</span>Here is the script as it is
now:</p>
<p>#!/usr/bin/python</p>
<p>from osgeo import gdal, osr,
gdalconst</p>
<p>import os, sys, time</p>
<p>import struct</p>
<p>import numpy as np</p>
<p>np.set_printoptions(threshold='nan')</p>
<p>gdal.AllRegister()</p>
<p>print 'Raster Source
1------------------'</p>
<p>ds1 =
gdal.Open('C:\Data\TE300by300.img',
gdal.GA_ReadOnly)</p>
<p>cols1 = ds1.RasterXSize</p>
<p>rows1 = ds1.RasterYSize</p>
<p>bands1 = ds1.RasterCount</p>
<p>print "Columns: " + str(cols1)</p>
<p>print "Rows: " + str(rows1)</p>
<p>print "Bands: " + str(bands1)</p>
<p>gt1 = ds1.GetGeoTransform()</p>
<p>width = ds1.RasterXSize</p>
<p>height = ds1.RasterYSize</p>
<p>minx = gt1[0]</p>
<p>print "Left(minx): "+ str(minx)</p>
<p>miny = gt1[3] + width*gt1[4] +
height*gt1[5] </p>
<p>print "Bottom(miny): "+ str(miny)</p>
<p>maxx = gt1[0] + width*gt1[1] +
height*gt1[2]</p>
<p>print "Right(maxx): "+ str(maxx)</p>
<p>maxy = gt1[3] </p>
<p>print "Top(maxy): "+ str(maxy)</p>
<p>pixWidth = gt1[1]</p>
<p>print "Pixel Width " +
str(pixWidth)</p>
<p>pixHeight = gt1[5]</p>
<p>print "Pixel Height " +
str(pixHeight)</p>
<p>xOrigin = gt1[0]</p>
<p>yOrigin = gt1[3]</p>
<p>print 'Raster Source
2------------------'</p>
<p>ds2 =
gdal.Open('C:\Data\LowElev300by300.img',
gdal.GA_ReadOnly)</p>
<p>cols2 = ds2.RasterXSize</p>
<p>rows2 = ds2.RasterYSize</p>
<p>bands2 = ds2.RasterCount</p>
<p>print "Columns: " + str(cols2)</p>
<p>print "Rows: " + str(rows2)</p>
<p>print "Bands: " + str(bands2)</p>
<p>gt2 = ds2.GetGeoTransform()</p>
<p>width = ds2.RasterXSize</p>
<p>height = ds2.RasterYSize</p>
<p>minx = gt2[0]</p>
<p>print "Left(minx): "+ str(minx)</p>
<p>miny = gt2[3] + width*gt2[4] +
height*gt2[5] </p>
<p>print "Bottom(miny): "+ str(miny)</p>
<p>maxx = gt2[0] + width*gt2[1] +
height*gt2[2]</p>
<p>print "Right(maxx): "+ str(maxx)</p>
<p>maxy = gt2[3] </p>
<p>print "Top(maxy): "+ str(maxy)</p>
<p>pixWidth = gt2[1]</p>
<p>print "Pixel Width " +
str(pixWidth)</p>
<p>pixHeight = gt2[5]</p>
<p>print "Pixel Height " +
str(pixHeight)</p>
<p>xOrigin = gt2[0]</p>
<p>yOrigin = gt2[3]</p>
<p> </p>
<p> </p>
<p>format = "HFA"</p>
<p>dst_file = "out.img"</p>
<p>dst_driver =
gdal.GetDriverByName(format)</p>
<p>dst_ds =
dst_driver.Create(dst_file, width,
height, 1, gdal.GDT_Byte )
#driver.Create( outfile, outwidth,
outheight, numbands, gdaldatatype)</p>
<p>empty = np.empty([height, width],
np.uint8 ) </p>
<p>srs = osr.SpatialReference()</p>
<p>dst_ds.SetProjection(ds2.GetProjection())</p>
<p>dst_ds.SetGeoTransform(ds2.GetGeoTransform())</p>
<p>dst_ds.GetRasterBand(1).WriteArray(empty)</p>
<p> </p>
<p> </p>
<p>#This is a way to go through a
given raster band </p>
<p>#one-by-one as an array by row,
getting all of the columns for</p>
<p>for r in range(rows1):</p>
<p><span> </span>data1 =
ds1.GetRasterBand(1).ReadAsArray(0,
r, cols1, 1)</p>
<p><span> </span>data1[data1>0]=1</p>
<p><span> </span>print "data1: "
+ str(data1)</p>
<p><span> </span>data2 =
ds2.GetRasterBand(1).ReadAsArray(0,
r, cols2, 1)</p>
<p><span> </span>data2[data2>0]=1</p>
<p><span> </span>print "data2: "
+ str(data2)</p>
<p><span> </span>result_bools =
np.logical_and(np.logical_and(data1
!= 0, data2 != 0), data1 == data2)<span>
</span></p>
<p><span> </span>result_ints =
np.array(result_bools, dtype=int)</p>
<p><span> </span>print
"result_ints: " + str(result_ints)</p>
<p><span> </span>dst_ds.GetRasterBand(1).WriteArray(result_ints,
0,
r)</p>
<p>dst_ds = None</p>
<p class="MsoNormal"> </p>
<p class="MsoNormal">Cheers, </p>
<p class="MsoNormal">Derek</p>
</div>
<br>
_______________________________________________<br>
gdal-dev mailing list<br>
<a moz-do-not-send="true"
href="mailto:gdal-dev@lists.osgeo.org"
target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a moz-do-not-send="true"
href="http://lists.osgeo.org/mailman/listinfo/gdal-dev"
target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote>
</div>
<br>
<br clear="all">
<br>
-- <br>
Best regards,<br>
Chaitanya kumar CH.<br>
<br>
<a moz-do-not-send="true"
href="tel:%2B91-9494447584"
value="+919494447584" target="_blank">+91-9494447584</a><br>
17.2416N 80.1426E<br>
</blockquote>
<br>
<br>
</div>
</div>
<span><font color="#888888">
<pre cols="72">--
Derek @ NEMAC</pre>
</font></span></div>
</blockquote>
</div>
<br>
<br clear="all">
<br>
-- <br>
Best regards,<br>
Chaitanya kumar CH.<br>
<br>
<a moz-do-not-send="true" href="tel:%2B91-9494447584"
value="+919494447584" target="_blank">+91-9494447584</a><br>
17.2416N 80.1426E<br>
</blockquote>
<br>
<br>
</div>
</div>
<span class="HOEnZb"><font color="#888888">
<pre cols="72">--
Derek @ NEMAC</pre>
</font></span></div>
</blockquote>
</div>
<br>
<br clear="all">
<br>
-- <br>
Best regards,<br>
Chaitanya kumar CH.<br>
<br>
+91-9494447584<br>
17.2416N 80.1426E<br>
</blockquote>
<br>
<br>
<pre class="moz-signature" cols="72">--
Derek @ NEMAC</pre>
</body>
</html>