3

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)
Improve this question
6
  • 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. Commented May 19, 2018 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. Commented May 19, 2018 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? Commented May 19, 2018 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. Commented May 19, 2018 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. Commented May 19, 2018 at 21:13

1 Answer 1

Reset to default
2

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.

Improve this answer
3
  • 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? Commented May 21, 2018 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. Commented May 21, 2018 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. Commented May 21, 2018 at 22:53

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.