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 &quot;data1: &quot; + 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 &quot;data2: &quot; + str(data2)</span><br style="font-family:courier new,monospace"><span style="font-family:courier new,monospace">    result_bools = np.logical_and((data1 &gt; 0), (data2 &gt; 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 &quot;result_ints: &quot; + 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">&lt;<a href="mailto:jdmorgan@unca.edu">jdmorgan@unca.edu</a>&gt;</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&gt;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&#39;t expect.<br>
    <br>
    Thanks, <br>
    Derek<div><div class="h5"><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&gt;</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&gt;</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">&lt;<a href="mailto:jdmorgan@unca.edu" target="_blank">jdmorgan@unca.edu</a>&gt;</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=&#39;nan&#39;)</p>
            <p>gdal.AllRegister()</p>
            <p>print &#39;Raster Source 1------------------&#39;</p>
            <p>ds1 = gdal.Open(&#39;C:\Data\TE300by300.img&#39;,
              gdal.GA_ReadOnly)</p>
            <p>cols1 = ds1.RasterXSize</p>
            <p>rows1 = ds1.RasterYSize</p>
            <p>bands1 = ds1.RasterCount</p>
            <p>print &quot;Columns: &quot; + str(cols1)</p>
            <p>print &quot;Rows: &quot; + str(rows1)</p>
            <p>print &quot;Bands: &quot; + 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 &quot;Left(minx): &quot;+ str(minx)</p>
            <p>miny = gt1[3] + width*gt1[4] + height*gt1[5] </p>
            <p>print &quot;Bottom(miny): &quot;+ str(miny)</p>
            <p>maxx = gt1[0] + width*gt1[1] + height*gt1[2]</p>
            <p>print &quot;Right(maxx): &quot;+ str(maxx)</p>
            <p>maxy = gt1[3] </p>
            <p>print &quot;Top(maxy): &quot;+ str(maxy)</p>
            <p>pixWidth = gt1[1]</p>
            <p>print &quot;Pixel Width &quot; + str(pixWidth)</p>
            <p>pixHeight = gt1[5]</p>
            <p>print &quot;Pixel Height &quot; + str(pixHeight)</p>
            <p>xOrigin = gt1[0]</p>
            <p>yOrigin = gt1[3]</p>
            <p>print &#39;Raster Source 2------------------&#39;</p>
            <p>ds2 = gdal.Open(&#39;C:\Data\LowElev300by300.img&#39;,
              gdal.GA_ReadOnly)</p>
            <p>cols2 = ds2.RasterXSize</p>
            <p>rows2 = ds2.RasterYSize</p>
            <p>bands2 = ds2.RasterCount</p>
            <p>print &quot;Columns: &quot; + str(cols2)</p>
            <p>print &quot;Rows: &quot; + str(rows2)</p>
            <p>print &quot;Bands: &quot; + 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 &quot;Left(minx): &quot;+ str(minx)</p>
            <p>miny = gt2[3] + width*gt2[4] + height*gt2[5] </p>
            <p>print &quot;Bottom(miny): &quot;+ str(miny)</p>
            <p>maxx = gt2[0] + width*gt2[1] + height*gt2[2]</p>
            <p>print &quot;Right(maxx): &quot;+ str(maxx)</p>
            <p>maxy = gt2[3] </p>
            <p>print &quot;Top(maxy): &quot;+ str(maxy)</p>
            <p>pixWidth = gt2[1]</p>
            <p>print &quot;Pixel Width &quot; + str(pixWidth)</p>
            <p>pixHeight = gt2[5]</p>
            <p>print &quot;Pixel Height &quot; + str(pixHeight)</p>
            <p>xOrigin = gt2[0]</p>
            <p>yOrigin = gt2[3]</p>
            <p> </p>
            <p> </p>
            <p>format = &quot;HFA&quot;</p>
            <p>dst_file = &quot;out.img&quot;</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&gt;0]=1</p>
            <p><span>    </span>print
              &quot;data1: &quot; + str(data1)</p>
            <p><span>    </span>data2 =
              ds2.GetRasterBand(1).ReadAsArray(0, r, cols2, 1)</p>
            <p><span>    </span>data2[data2&gt;0]=1</p>
            <p><span>    </span>print
              &quot;data2: &quot; + 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
              &quot;result_ints: &quot; + 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 href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
          <a 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 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>