7

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.

3

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.

  • This fails at the file.GetRasterBand(1).WriteArray(result) line ("AttributeError: 'NoneType' object has no attribute 'GetRasterBand'"). Also, could you comment the code to explain how it works? – Simbamangu Aug 30 '16 at 15:43
  • @Simbamangu It's strange, because I've laready generated raster with this code, maybe something is wrong with your QGIS Python installation. – dmh126 Aug 30 '16 at 15:46
  • @dmh126 tried it and got the sam error message as Simbamangu – Kurt Aug 30 '16 at 16:52
  • @Kurt OK I agree, same problem on Windows, but on Linux it works good. I'll try to figure this out. – dmh126 Aug 30 '16 at 17:00
  • OK now it should works - just add gdal.UseExceptions() and edit output_path to path where you have permissions to write files. I've done this on Windows 10 and it works. – dmh126 Aug 30 '16 at 17:14
1

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()

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.