Creating raster where each cell records distance to sea?

Question

I want to create a raster with a 25 metre × 25 metre resolution, where each cell contains the distance to the nearest coastline, as calculated from the center of the cell. To do this, all I have is a shapefile of the coastlines of New Zealand.

I have tried following Dominic Roye's tutorial for doing it in R which works... kind of. It is fine down to about a 1 km × 1 km resolution but if I try to go any higher the RAM it requires well exceeds that available on my PC (~ 70 gb of RAM required) or any other that I have access too. In saying that, I think this is a limitation of R and I suspect that QGIS might have a more computationally efficient way of creating this raster, but I am new to it and I can't quite figure out how to do it.

I have tried following Creating raster with distance to feature using QGIS? to create it in QGIS but it returns this error:

_core.QgsProcessingException: Could not load source layer for INPUT: C:/..../Coastline/nz-coastlines-and-islands-polygons-topo-150k.shp not found

and I am not sure why.

Does anyone have any suggestions of what might be going wrong or an alternative way of doing this?

Edit:

The raster I am hoping to produce would have about 59684 rows and 40827 columns so that it overlaps with the annual water deficit raster from LINZ. If the raster that is produced is larger than the annual water deficit raster I can snip it in R though...

One thing that I think might be a potential issue is that the shapefile of NZ's coastline has a large amount of sea between the islands, and I am not interested in calculating the distance to coast for these cells. I really only want to calculate the values for cells that are include some slice of land. I am not sure how to do this, or if it is an actually problem though.

Answer 1

With PyQGIS and GDAL python library is not very difficult to do that. You need geo transform parameters (top left x, x pixel resolution, rotation, top left y, rotation, n-s pixel resolution) and rows and columns number for creating resulting raster. For calculating the distance to the nearest coastline, it is necessary a vector layer for representing coastline.

With PyQGIS, each raster point as center of the cell is calculated and its distance to coastline is measured by using 'closestSegmentWithContext' method from QgsGeometry class. GDAL python library is used for producing a raster with these distance values in rows x columns array.

Following code was used for creating a distance raster (25 m × 25 m resolution and 1000 rows x 1000 columns) starting in point (397106.7689872353, 4499634.06675821); near to west coastline of USA.

from osgeo import gdal, osr
import numpy as np
from math import sqrt

registry = QgsProject.instance()

line = registry.mapLayersByName('shoreline_10N')

crs = line[0].crs()

wkt = crs.toWkt()

feats_line = [ feat for feat in line[0].getFeatures()]

pt = QgsPoint(397106.7689872353, 4499634.06675821)

xSize = 25
ySize = 25

rows = 1000
cols = 1000

raster = [ [] for i in range(cols) ]

x =   xSize/2
y = - ySize/2

for i in range(rows):
    for j in range(cols):
        point = QgsPointXY(pt.x() + x, pt.y() + y)
        tupla = feats_line[0].geometry().closestSegmentWithContext(point)
        raster[i].append(sqrt(tupla[0]))

        x += xSize
    x =  xSize/2
    y -= ySize

data = np.array(raster)

# Create gtif file 
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/distance_raster.tif"

dst_ds = driver.Create(output_file, 
                       cols, 
                       rows, 
                       1, 
                       gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( data )

transform = (pt.x(), xSize, 0, pt.y(), 0, -ySize)

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(transform)

# setting spatial reference of output raster 
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )

dst_ds = None

After running above code, resulting raster was loaded in QGIS and it looks as in following image (pseudocolor with 5 classes and Spectral ramp). Projection is UTM 10 N (EPSG:32610)

Answer 2

May be a solution to try:

1. Generate a grid (Type "point", algorithm "Create grid")
2. Calculate the nearest distance between your points (grid) and your line (coast) with the algorithm "join attribute by nearest". Be carefull to choose only a maximum of 1 nearest neighbors.

Now you should have a new point layer with the distance to the coast like in this example

3. If needed, you can convert your new point layer to a raster (algorithm "rasterize")

Answer 3

Within QGIS you could try the GRASS plugin. As far as I know it manages the memory better than R, and I expect the other solution to fail on large areas.

the GRASS command is called r.grow.distance , which you can find in the processing toolbar. Note that you need to convert your line to raster at first.

One of your issue could be the size of the output, so you can add some usefull creation options such as (for a tif file) BIGTIFF=YES,TILED=YES,COMPRESS=LZW,PREDICTOR=3

Answer 4

I wasn't originally going to answer my own question but a colleague of mine (who does not use this site), wrote me a bunch of python code to do what I am after; including limiting the cells to have the distance to coast for only terrestrial cells and leaving the sea based cells as NAs. The following code should be able to run from any python console, with the only things needing alteration being:

1) Put the script file in the same folder as the shape file of interest;

2) change the name of the shapefile in the python script to whatever the name of your shapefile is;

3) set the desired resolution, and;

4) change the extent to match other rasters.

Larger shapefiles than what I am using will require large amounts of RAM but otherwise the script is quick to run (about three minutes to produce a 50m resolution raster and ten minutes for 25m resolution raster).

#------------------------------------------------------------------------------

from osgeo import gdal, ogr
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
import time

startTime = time.perf_counter()

#------------------------------------------------------------------------------

# Define spatial footprint for new raster
cellSize = 50 # ANDRE CHANGE THIS!!
noData = -9999
xMin, xMax, yMin, yMax = [1089000, 2092000, 4747000, 6224000]
nCol = int((xMax - xMin) / cellSize)
nRow = int((yMax - yMin) / cellSize)
gdal.AllRegister()
rasterDriver = gdal.GetDriverByName('GTiff')
NZTM = 'PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],AUTHORITY["EPSG","2193"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

#------------------------------------------------------------------------------ 

inFile = "new_zealand.shp" # CHANGE THIS!!

# Import vector file and extract information
vectorData = ogr.Open(inFile)
vectorLayer = vectorData.GetLayer()
vectorSRS = vectorLayer.GetSpatialRef()
x_min, x_max, y_min, y_max = vectorLayer.GetExtent()

# Create raster file and write information
rasterFile = 'nz.tif'
rasterData = rasterDriver.Create(rasterFile, nCol, nRow, 1, gdal.GDT_Int32, options=['COMPRESS=LZW'])
rasterData.SetGeoTransform((xMin, cellSize, 0, yMax, 0, -cellSize))
rasterData.SetProjection(vectorSRS.ExportToWkt())
band = rasterData.GetRasterBand(1)
band.WriteArray(np.zeros((nRow, nCol)))
band.SetNoDataValue(noData)
gdal.RasterizeLayer(rasterData, [1], vectorLayer, burn_values=[1])
array = band.ReadAsArray()
del(rasterData)

#------------------------------------------------------------------------------

distance = ndimage.distance_transform_edt(array)
distance = distance * cellSize
np.place(distance, array==0, noData)

# Create raster file and write information
rasterFile = 'nz-coast-distance.tif'
rasterData = rasterDriver.Create(rasterFile, nCol, nRow, 1, gdal.GDT_Float32, options=['COMPRESS=LZW'])
rasterData.SetGeoTransform((xMin, cellSize, 0, yMax, 0, -cellSize))
rasterData.SetProjection(vectorSRS.ExportToWkt())
band = rasterData.GetRasterBand(1)
band.WriteArray(distance)
band.SetNoDataValue(noData)
del(rasterData)

#------------------------------------------------------------------------------

endTime = time.perf_counter()

processTime = endTime - startTime

print(processTime)

Answer 5

I would try other way around. If you are using poligon of NZ then convert polygon edges to line. After that create buffer on the boundary for every 25 meters of distance from boundary(maybe centorid might help in determing when to stop). Then cut buffers out with polygon and then convert that polygons to raster. I am not sure this would work but definitely you gonna need less RAM. And PostGiS is great when you have performance issues.

Hope it might help at least a little bit :)