10

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.

  • 1
    Are you running a script to do this? Or are you using the tools in QGIS? Something to check, even though it sounds like it should - check the file actually exists where you say it does...also check that you have read and write access to that particular folder. – Keagan Allan Sep 19 at 23:39
  • Currently using the tools but I am quite keen to learn the script, just not sure where to start. I am sure the file exists, as I have loaded the .shp file into QGIS and it pops up as an image. I should have read/write access too as I am an admin on the machine and it is just in my dropbox. – André.B Sep 19 at 23:42
  • Try moving it out of Dropbox to a local drive. There may be an issue with the path causing QGIS to reject it. What you are looking to do should be pretty simple in QGIS. Which version of QGIS are you using? – Keagan Allan Sep 19 at 23:45
  • 1
    Ok, try converting the polyline to a raster. The Proximity tool in QGIS needs a raster input. Play around with the settings as per the tool's help: docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/…. Take note, it is still an intensive process, I am testing it for fun now and it has been running for 30mins and still going... – Keagan Allan Sep 20 at 0:29
  • 1
    What size of output raster in terms of rows and columns are you trying to create? Are you actually going to be able to work with that raster once you do create it? If the file size of the whole thing is a problem, could you create smaller tiles, which is also something you can do in parallel on a cluster or cloud for speed. – Spacedman Sep 22 at 16:49
9
+50

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)

enter image description here

  • This might not be an issue but one thing that I am a little worried about is that the polygon is of New Zealand and its surrounding islands which means it includes a huge amount of the surrounding sea. I am trying to get my head around the code, but with your example would I be able to set the value for all the cells at sea to NA? I am really only interested in the distance to the sea from points on land. – André.B Sep 20 at 2:56
  • Apologies in advance if this is a dumb question but how do I choose a new starting point in New Zealand in the way that you set the coordinates for the states? Also how do I keep it in EPSG:2193? – André.B Sep 23 at 1:41
6

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 enter image description here

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

enter image description here

2

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.

enter image description here

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

  • Is there a way that I could eliminate any sea area so as to reduce the size/computation time? – André.B Sep 25 at 4:20
  • in theory, if you use the distance from the see (all see pixels with same value, that is a polygon) instead of from the coast line you should save computation time. The uncompressed size of the raster should be the same, but compression will be more efficient, so you should also reduce the final size. – radouxju Sep 25 at 6:13
0

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

0

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)

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.