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> 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> 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 >> http://lists.osgeo.org/mailman/listinfo/gdal-dev >> > > > > -- > Best regards, > Chaitanya kumar CH. > > +91-9494447584 > 17.2416N 80.1426E > > > > -- > Derek @ NEMAC > > -- Best regards, Chaitanya kumar CH. +91-9494447584 17.2416N 80.1426E
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev