2

As I work through learning programing with raster imagery, I am trying to read in a raster, change the pixel values and write the output file. Below is the code. I have blocked out the change pixel values as my output file was empty and in the wrong location, so I am now trying to just read in a file by block and write out the same file. However the file I get is defaulted to the 64/64 block size, and my skip factor is 999, so my output raster does not 'match' my input raster. I have not been able to SET the metadata changes. Can someone tell me the correct way in a script to set new metadata options, and then I will need the numpy.where to process correctly?

import os, sys, ogr, gdal, utils, numpy
from gdalconst import *

file = (r'C:\Users\Desktop\images\large_raster.img')

gdal.AllRegister()

ds = gdal.Open(file, GA_ReadOnly)
if ds is None:
    print 'Could not open file'
    sys.exit(1)

rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

band = ds.GetRasterBand(1)

yBlocksize = 512
xBlocksize = 512
count = 0

driver = gdal.GetDriverByName('HFA')
metadata = driver.GetMetadata()

#print metadata

outfile = (r'C:\Users\aschilli\Desktop\images\chrisG_test2.img')

outDataset = driver.Create(outfile, cols, rows, 1)

outBand = outDataset.GetRasterBand(1)
outBand.FlushCache()
outBand.SetNoDataValue(-99) #this gives 0
outDataset.SetGeoTransform(ds.GetGeoTransform())
outDataset.SetProjection(ds.GetProjection()) 
gdalinfo = (r"C:\Python27\ArcGIS10.3\Scripts\new_ve_folder\Lib\site-packages\osgeo\gdalinfo.exe")
command = gdalinfo + ' -hist ' + ' -stats ' + outfile
os.system(command)
gdaladdo = (r"C:\Python27\ArcGIS10.3\Scripts\new_ve_folder\Lib\site-packages\osgeo\gdaladdo.exe")
command2 = gdaladdo + ' --config ' ' HFA_USE_RRD ' ' YES ' + outfile + ' 3 9 27 81'
os.system(command2)
#this code below does nothing
outDataset.SetMetadataItem({
    'BLOCKSIZE': '512/512',
    'STATISTICS': 'TRUE',
    'DMD_CREATIONDATATYPES': 'Int16'
})

#loop through rows
for i in range(0, rows, yBlocksize):
    if i + yBlocksize < rows:
        numRows = yBlocksize
    else:
        numRows = rows - i
#loop through coloumns
    for j in range(0, cols, xBlocksize):
        if j + xBlocksize < cols:
            numCols = xBlocksize
        else:
            numCols = cols - j

    #read data and do calculations
        data = band.ReadAsArray(j, i,  numCols, numRows)
        # out_data = numpy.where(abs(data) < 500, 10, data)

        outBand.WriteArray(data, j, i)

del outDataset

#http://gis.stackexchange.com/questions/160578/write-python-ndarray-to-raster?'
1

I have figured out and gotten to work one part of my statistics issue. src_ds = gdal.Open(outfile) outBand.SetMetadataItem("LAYER_TYPE", "thematic") outBand.SetMetadataItem("STATISTICS", "YES")

This will give me my skip factor of 1 for ERDAS img files.

My file only writes out the first block of data, so I am missing something telling my script to write out BY block and not just the first block.

  • adding < outBand.SetMetadataItem("BLOCKSIZE", "512/512")> solved the problem of only writing first block. I have gotten most of this to work, now for some really detailed reclassifying. – als May 12 '17 at 15:41
0

Your question is a bit convoluted. I gather that you are trying to classify an image using a numpy.where statement, and are having problems setting the correct georeferencing on the output image. From your script, I can also see that you are reading only the first raster band in GDAL, so I am assuming that your input raster has only one band.

The following code will accomplish what you want. The commands OutRaster.SetProjection() and OutRaster.SetGeoTransform() set the georeferencing of your output raster:

from osgeo import gdal
import numpy

InputImage = 'ImageName.tif'
OutImage = 'OutImageName.tif'
gdalformat = 'GTiff'
########################################################################
print("Opening "+InputImage+"...")
Dataset = gdal.Open(InputImage, gdal.GA_ReadOnly)
XSize = Dataset.RasterXSize
YSize = Dataset.RasterYSize
Projection = Dataset.GetProjectionRef()
GeoTransform = Dataset.GetGeoTransform()

Band1 = Dataset.GetRasterBand(1) # Select band 1 to read as numpy array
Array = Band1.ReadAsArray()
Dataset = None # Closes the input image

# Create an empty duplicate numpy array.  
Classification = numpy.empty_like(Array, dtype='uint8')

print("Perfoming classification...")
Classification = numpy.where((Array <= 500), 1, Classification)
Classification = numpy.where((Classification != 1) & (Array > 500), 2, Classification)
Array = None # Close the input array

# Create the output raster
driver = gdal.GetDriverByName(gdalformat)
metadata = driver.GetMetadata()
OutRaster = driver.Create(OutImage, XSize, YSize, 1, gdal.GDT_Byte) # (OutImageName, XSize, YSize, Image bands, datatype)
OutRaster.SetProjection(Projection)
OutRaster.SetGeoTransform(GeoTransform)

# Write classification to band 1
OutBand = OutRaster.GetRasterBand(1)
OutBand.SetNoDataValue(0)
OutBand.WriteArray(Classification)

# Close the datasets
OutRaster, OutBand, Classification = None, None, None
print("Done.")

I don't know why you are trying to read the raster in blocks - perhaps because the raster is very large? In this case, you could tile the image and classify each tile before mosaicking them back together. Alternatively, you could use something like rios, which can read the raster in chunks rather than loading the entire image into memory.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.