Derek, You can assign a color table. Ex: (0,0,0) for 0, (255,255,255) for 1. You can assign a color table to a raster band using the method SetRasterColorTable(). You can find some sample python code in GDAL's test suite[1][2].
[1]: http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/colortable.py [2]: http://trac.osgeo.org/gdal/browser/trunk/autotest/gcore/tiff_write.py#L916 On Thu, Feb 23, 2012 at 7:45 PM, jdmorgan <jdmor...@unca.edu> wrote: > Thanks Chaitanya, > > Actually this seems to be working, which is pretty cool. I think where I > am getting confused is that I am attempting to verify the results by > loading the resulting raster data into a GIS e.g. ArcGIS or QGIS and doing > spot checking. However, as I haven’t done any actual color coding on the > classes e.g. making the 1 black and 0’s white, I think I am not able to > verify it correctly. If I let ArcGIS do a default classification it does > something line 0’s white and 0-1 black. Is there a way to color class > this in GDAL python? That would probably be a better approach. > > Cheers, > > Derek > > > On 2/22/2012 10:01 AM, Chaitanya kumar CH wrote: > > 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 > > > > -- > 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
