I am creating grayscale, single channel geotiff files from numpy arrays with the below module, but I would like to apply a fairly standard colormap to the data. I've tried three channels of R,B and G and have had no success if producing anything that looks correct. I'm thinking my approach, in general is flawed. Any points in the right direction would be appreciated. How can I go from grayscale to color geotiffs?
def geotiff_output(outfile, mag_grid, lons, lats): #*********************Test*************************** from osgeo import gdal, osr, ogr xres = lons[1] - lons[0] yres = lats[1] - lats[0] ysize = len(lats) xsize = len(lons) ulx = lons[0] - (xres / 2.) uly = lats[-1] - (yres / 2.) driver = gdal.GetDriverByName('GTiff') ds = driver.Create(sc_settings.dynamic_folder_script + outfile[:-4] + ".tif", xsize, ysize, 1, gdal.GDT_Byte, ) # this assumes the projection is Geographic lat/lon WGS 84 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) ds.SetProjection(srs.ExportToWkt()) gt = [ulx, xres, 0, uly, 0, yres ] ds.SetGeoTransform(gt) outband=ds.GetRasterBand(1) #mask no data points... mag_grid1 = np.ma.masked_where(mag_grid==-99.,mag_grid) #assign actual data values for inclusion in metadata outband.SetStatistics(np.round(np.min(mag_grid1),decimals=4), np.round(np.max(mag_grid1),decimals=4), np.round(np.mean(mag_grid1),decimals=4), np.round(np.std(mag_grid1),decimals=4)) #normalize the array to a 0-255 scale if np.min(mag_grid1) < 0: mag_grid1 = mag_grid1 + np.abs(np.min(mag_grid1)) #Normalize array to a 0-255 scale for raster channel mag_grid2 = ((mag_grid1 - np.min(mag_grid1))/(np.max(mag_grid1) - np.min(mag_grid1)))*256 outband.WriteArray(mag_grid2) ds = None
أكثر...
def geotiff_output(outfile, mag_grid, lons, lats): #*********************Test*************************** from osgeo import gdal, osr, ogr xres = lons[1] - lons[0] yres = lats[1] - lats[0] ysize = len(lats) xsize = len(lons) ulx = lons[0] - (xres / 2.) uly = lats[-1] - (yres / 2.) driver = gdal.GetDriverByName('GTiff') ds = driver.Create(sc_settings.dynamic_folder_script + outfile[:-4] + ".tif", xsize, ysize, 1, gdal.GDT_Byte, ) # this assumes the projection is Geographic lat/lon WGS 84 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) ds.SetProjection(srs.ExportToWkt()) gt = [ulx, xres, 0, uly, 0, yres ] ds.SetGeoTransform(gt) outband=ds.GetRasterBand(1) #mask no data points... mag_grid1 = np.ma.masked_where(mag_grid==-99.,mag_grid) #assign actual data values for inclusion in metadata outband.SetStatistics(np.round(np.min(mag_grid1),decimals=4), np.round(np.max(mag_grid1),decimals=4), np.round(np.mean(mag_grid1),decimals=4), np.round(np.std(mag_grid1),decimals=4)) #normalize the array to a 0-255 scale if np.min(mag_grid1) < 0: mag_grid1 = mag_grid1 + np.abs(np.min(mag_grid1)) #Normalize array to a 0-255 scale for raster channel mag_grid2 = ((mag_grid1 - np.min(mag_grid1))/(np.max(mag_grid1) - np.min(mag_grid1)))*256 outband.WriteArray(mag_grid2) ds = None
أكثر...