How can I generate linear elevation raster?

Question

I want to generate a raster for testing purposes. I'm using QGIS and I was hoping it would have the tools to do this but I couldn't find a way.

The raster should be 1001 m x 1001 m in size with each point being 1 m x 1 m. The center point should be 500m high and the surrounding points (including the diagonal points) should linearly descent until reaching 0m at the very edge of the pyramid (i.e. 1m descent per point).

This raster should be centered to geographic location 0.0, 0.0 in WGS84 datum.

I'm hoping there's a way to do this in QGIS but other (free) tools are welcome as well.

Answer 1

With Python script (just paste it into the Python console in QGIS, or run as a .py file there):

import numpy as np
from osgeo import gdal, osr
gdal.UseExceptions()

x = np.arange(0,500,1)
y = np.arange(0,500,1)
xx,yy = np.meshgrid(x,y)
l = np.tril(xx)
u = np.triu(yy,k=1)
m1 = l+ u #lu
m2 = np.rot90(m1,1) #ll
m3 = np.rot90(m1,2) #rl
m4 = np.rot90(m1,3) #ru

h1=np.concatenate((m1,m4), axis=1)
h2=np.concatenate((m2,m3), axis=1)

result = np.concatenate((h1,h2))

driver = gdal.GetDriverByName('GTiff')
output_path = "D:\\raster.tif"
file = driver.Create(output_path, 1000 ,1000 , 1, gdal.GDT_Float32)
file.GetRasterBand(1).WriteArray(result)

# spatial reference goes here
proj = osr.SpatialReference()
proj.ImportFromEPSG(4326)   
file.SetProjection(proj)
file.SetGeoTransform((0.0, 0.005, 0.0, 0.0, 0.0, -0.005))

file.FlushCache()

Some mathematics should be done. I hope it works.

Answer 2

Solution given by dmh126 is close but not quite what I needed (the center point is not at 0, 0 and the raster is not 1001m x 1001m and the center height is not 500m). I modified his answer a bit and here is what I did:

import numpy as np
from osgeo import gdal
gdal.UseExceptions()

x = np.linspace(0,500,500)
y = np.linspace(0,500,500)
xx,yy = np.meshgrid(x,y)
l = np.tril(xx)
u = np.triu(yy,k=1)
m1 = l+ u #lu
m2 = np.rot90(m1,1) #ll
m3 = np.rot90(m1,2) #rl
m4 = np.rot90(m1,3) #ru

h1=np.concatenate((m1,m4), axis=1)
h2=np.concatenate((m2,m3), axis=1)

result = np.concatenate((h1,h2))

driver = gdal.GetDriverByName('GTiff')
output_path = "C:\\temp\\raster.tif"
file = driver.Create(output_path, 1001 ,1001 , 1, gdal.GDT_Float32)
gt = file.GetGeoTransform()
gtl = list(gt)
gtl[0] -= 500.0
gtl[3] -= 500.0
file.SetGeoTransform(tuple(gtl))
file.GetRasterBand(1).WriteArray(result)
file.FlushCache()