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 <jdmor...@unca.edu
<mailto:jdmor...@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 <jdmor...@unca.edu
<mailto:jdmor...@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@lists.osgeo.org <mailto:gdal-dev@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
