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" class="br0">[</span><span style="font-family:courier new,monospace">data1></span><span style="font-family:courier new,monospace" class="nu0">0</span><span style="font-family:courier new,monospace" class="br0">]</span><span style="font-family:courier new,monospace">=</span><span style="font-family:courier new,monospace" class="nu0">1<br>
...<br></span><span style="font-family:courier new,monospace">data2</span><span style="font-family:courier new,monospace" class="br0">[</span><span style="font-family:courier new,monospace">data2></span><span style="font-family:courier new,monospace" class="nu0">0</span><span style="font-family:courier new,monospace" class="br0">]</span><span style="font-family:courier new,monospace">=</span><span style="font-family:courier new,monospace" class="nu0">1</span><br>
<br><br><div class="gmail_quote">On Wed, Feb 22, 2012 at 8:46 AM, jdmorgan <span dir="ltr"><<a 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">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><u></u> <u></u></p>
<p><u></u> <u></u></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><u></u> <u></u></p>
<p><u></u> <u></u></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"><u></u> <u></u></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">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>+91-9494447584<br>17.2416N 80.1426E<br>