Mdthbs

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?

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)

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