1

What I want to realize is that

1, read a raster file as numpy array;

2, perform Gaussian blur for the array;

3, convert the array to raster file.

When I convert the array to the raster file, I need the select4.tif to be the same Xmin, Xmax, Ymin, Ymax as select4.tif file. In which the function ReprojectCoords is to project the same coordinate as select3.tif to select4.tif, and GetExtent is to fetch the extent from the select3.tif so as to set the same extent value and raster size to select4.tif.

Here are respectively the select3, select4 raster files.

The original raster file looks like select3

After the processing and conversion, the select3.tif and select4.tif do not have the same extent enter image description here

I think the problem lies in the function GetExtent which the source codes are from this post How to get raster corner coordinates using Python GDAL bindings?. But I can't figure out why it generated a different extent as the original select3.tif file. Could anyone help me to find out where the problem is?

The following is the codes for the whole conversion procedure.

read the raster file as numpy array and do the Gaussian blur.

import numpy as np
import gdal, ogr, os, osr 
#gaussian blur
ds = gdal.Open("C:\\Users\\select\\select3.tif")
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())


# Gaussian blur

def gaussian_blur(in_array, size):
    # expand in_array to fit edge of kernel
    padded_array = np.pad(in_array, size, 'symmetric')
    # build kernel
    x, y = np.mgrid[-size:size + 1, -size:size + 1]
    g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
    g = (g / g.sum()).astype(in_array.dtype)
    # do the Gaussian blur
    return fftconvolve(padded_array, g, mode='valid')



Gblur= gaussian_blur(myarray, size=10) # return numpy array 

define a function to project the coordinate

def ReprojectCoords(coords,src_srs,tgt_srs):
    ''' Reproject a list of x,y coordinates.

        @type geom:     C{tuple/list}
        @param geom:    List of [[x,y],...[x,y]] coordinates
        @type src_srs:  C{osr.SpatialReference}
        @param src_srs: OSR SpatialReference object
        @type tgt_srs:  C{osr.SpatialReference}
        @param tgt_srs: OSR SpatialReference object
        @rtype:         C{tuple/list}
        @return:        List of transformed [[x,y],...[x,y]] coordinates
    '''
    trans_coords=[]
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    for x,y in coords:
        x,y,z = transform.TransformPoint(x,y)
        trans_coords.append([x,y])
    return trans_coords

to get the Xmin, Ymin of the raster file

def GetExtent(gt,cols,rows):
    ''' Return list of corner coordinates from a geotransform

        @type gt:   C{tuple/list}
        @param gt: geotransform
        @type cols:   C{int}
        @param cols: number of columns in the dataset
        @type rows:   C{int}
        @param rows: number of rows in the dataset
        @rtype:    C{[float,...,float]}
        @return:   coordinates of each corner
    '''
    ext=[]
    xarr=[0,cols]
    yarr=[0,rows]

    for px in xarr:
        for py in yarr:
            x=gt[0]+(px*gt[1])+(py*gt[2])
            y=gt[3]+(px*gt[4])+(py*gt[5])
            ext.append([x,y])
            #print (x,y)
        yarr.reverse()
    return ext
ds = gdal.Open("C:\\Users\\select\\select3.tif")

gt=ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
ext=GetExtent(gt,cols,rows)

project the spatial reference to the ext

src_srs = osr.SpatialReference()
src_srs.ImportFromWkt(ds.GetProjection())
tgt_srs=osr.SpatialReference()
tgt_srs.ImportFromEPSG(4326)
tgt_srs = src_srs.CloneGeogCS()
geo_ext = ReprojectCoords(ext, src_srs, tgt_srs)
print(geo_ext)
originXmin = geo_ext[0][0]

originYmin = geo_ext[1][1]

convert array to raster file

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):

    cols = array.shape[1]
    rows = array.shape[0]

    originX =rasterOrigin[0]  # 36.65808 
    originY =rasterOrigin[1]   # 44.36431

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()



def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    reversed_arr = array[::-1] # reverse array so the tif looks like the array
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster



array = Gblur
rasterOrigin = (originXmin,originYmin)
pixelWidth = 0.0001 
pixelHeight = 0.0001

newRasterfn = "C:\\Users\\select\\select4.tif"


main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)
  • add "outRaster = None" to the bottom of your array2raster function. I don't know if that's the solution but it might not be writing the new rasters to disk properly. I am also not sure why you are doing any reprojecting, as both your initial images already have the same CRS and extents. Load the array, manipulate it, write it back to disk using its original geotransform and CRS. No reprojection necessary. – Jon May 19 '18 at 20:32
  • Oh, and another thing that is probably causing problems--it looks like your Gaussian blur is returning an array of a different size than you input. This will make the original geotransform incorrect. The easiest way to handle it is to "unpad" your array before returning it from your Gaussian blur function. It's a one-liner using slicing. If you want to retain the padding, you'll have to adjust your x,y origin by subtracting xres * size, yres * size from each. – Jon May 19 '18 at 20:40
  • Thank you. What do you by ` you'll have to adjust your x,y origin by subtracting xressize, yressize from each`. In which place should I adjust the x, y origin? – statistics_learning May 19 '18 at 20:44
  • Your line gt=ds.GetGeoTransform() returns the geotransform of the original image, which is [127, 50]. But after you process your array with Gaussian blur, the returned array has been extended by size on each side. So your returned array is [147, 70] for size =10. This means that your upper-left coordinate has changed. The geotransform contains the x and y value of the original upper-left coordinate; you need to modify it before sending it to GetExtents to match your padded array. – Jon May 19 '18 at 20:53
  • Sorry, I am still not very clear since I call the GetExtents function before the Gaussian blur, when convert array after Gaussian blur to raster, I reset the coordinate based on the original image (which is the return values by GetExtents). So although the array may be modified by the Gaussian blur, I reset it again in the array2raster function. Could you write out a solution and highlight which line of codes should be modified if you already clear where the problem is. Thank you. – statistics_learning May 19 '18 at 21:13
1

I am not entirely sure where your problem is, but here is a fully-working solution you can use based on code I use frequently to write geotiffs. Note that this only works when you already have a geotiff whose georeferencing info you want to copy to the array (as you do--both s3 and s4 are geotiffs). When you call write_gtiff() here, just provide the gdal object for s3 when you are writing s4, as I've done in the code.

import numpy as np
import gdal, osr 
from scipy.signal import fftconvolve

def gaussian_blur(in_array, size):
    # expand in_array to fit edge of kernel
    padded_array = np.pad(in_array, size, 'symmetric')
    # build kernel
    x, y = np.mgrid[-size:size + 1, -size:size + 1]
    g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
    g = (g / g.sum()).astype(in_array.dtype)
    # do the Gaussian blur
    return fftconvolve(padded_array, g, mode='valid')


def write_gtiff(array, gdal_obj, outputpath, dtype=gdal.GDT_UInt16, options=0, color_table=0, nbands=1, nodata=False):
    """
    Writes a geotiff.

    array: numpy array to write as geotiff
    gdal_obj: object created by gdal.Open() using a tiff that has the SAME CRS, geotransform, and size as the array you're writing
    outputpath: path including filename.tiff
    dtype (OPTIONAL): datatype to save as
    nodata (default: FALSE): set to any value you want to use for nodata; if FALSE, nodata is not set
    """

    gt = gdal_obj.GetGeoTransform()

    width = np.shape(array)[1]
    height = np.shape(array)[0]

    # Prepare destination file
    driver = gdal.GetDriverByName("GTiff")
    if options != 0:
        dest = driver.Create(outputpath, width, height, nbands, dtype, options)
    else:
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:
        dest.GetRasterBand(1).SetColorTable(color_table)

    dest.GetRasterBand(1).WriteArray(array)

    if nodata is not False:
        dest.GetRasterBand(1).SetNoDataValue(nodata)

    # Set transform and projection
    dest.SetGeoTransform(gt)
    wkt = gdal_obj.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dest.SetProjection(srs.ExportToWkt())

    # Close output raster dataset 
    dest = None



s3 = r"C:\Users\Jon\Desktop\select3.tif"
s4 = r"C:\Users\Jon\Desktop\select4.tif"

ds3 = gdal.Open(s3)
ds4 = gdal.Open(s4)

array = np.array(ds3.GetRasterBand(1).ReadAsArray())
Gblur = gaussian_blur(array, size=10)
outpath = r"C:\Users\Jon\Desktop\test3.tif"
write_gtiff(Gblur, ds3, outpath)

array = np.array(ds4.GetRasterBand(1).ReadAsArray())
Gblur = gaussian_blur(array, size=10)
outpath = r"C:\Users\Jon\Desktop\test4.tif"
write_gtiff(Gblur, ds3, outpath)

The options variable is for creation specifications, like tile size, compression, etc. Stuff that you don't need to worry about for such a small geotiff.

I checked the result in QGIS and the extents are the same for test3, test4, and s3.

  • Thank you Jon, you did solved my problem. I go through your codes carefully but not very sure which lines of codes are the key factor that solves the problem. Cuz it seems there are a lot of changes. Could you indicate me how do you force the two raster files to give the same extents? – statistics_learning May 21 '18 at 21:49
  • The most important thing you need to understand is the GeoTransform (gt). It contains the coordinates of the upper-left-most pixel and the resolutions (x and y) of a geotiff. All I did was take the geotransform from s3 and assign it to the processed s4 (dest.SetGeoTransform(gt)). This only works because 3 and 4 have the same resolution and CRS. If you really want to understand, read more about gdal's geotransform property. – Jon May 21 '18 at 22:22
  • @statistics_learning Basically, there is no need for your ReprojectCoords or GetExtent functions. You can just copy the georeferencing information (CRS and geotransform) directly from s3 onto the output files. That's what my script is doing. In the future, if you have two geotiffs that have different CRS (i.e. EPSG) or different resolutions that you want to match, then you will need to worry about reprojecting. – Jon May 21 '18 at 22:53

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.