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
17.2416N 80.1426E
--
Derek @ NEMAC
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
