Creating raster where each cell records distance to sea?
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)
May be a solution to try:
- Generate a grid (Type "point", algorithm "Create grid")
- 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
- If needed, you can convert your new point layer to a raster (algorithm "rasterize")
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