Assign nearest raster cell value to point

Question

I came across to a similar issue described in Finding nearest raster cell value based on vector point?

I understand the proposed solution with Raster to Point and then Spatial Join with (for example) Closest Geodesic option. However, is not clear how to use the coordinate formula described by @whuber:

(x',y') = ( min(max(x,x0),x1), min(max(y,y0),y1 )

For example on my global raster x0=-180, x1=180, y0=-90 and y1=90 If I apply the above formula for each given (x,y) then I get no difference (no shift). That is obvious, because all the (x,y) points for which I want to extract value from raster are within the extent of the raster, just that they fall in NoData cells. I presume the formula makes sense when knowing the x0,x1,y0,y1 of the nearest cell that contains any valid value different from NoData (which I need to find in the first place).

Answer 1

you could create thiessen polygons around each raster value (using euclidian allocation tool). Then you can extract the value for each point using extract value to point.

Note that distance analysis at global scale has two problems: 1) there is a rupture in the map (usually between -180 and + 180) and the scale is not constant.

Answer 2

This is a more complex problem because the "coordinate formula" describe only an ideal case. If you have a raster clipped by a shapefile the formula will not work. However, the issue can be totally solve with Python code. To test my approach I used a little raster (20 x 20) to easily corroborate the results. For this situation, a shapefile was used to extract a mask with nodata values equals to -999.

I got the next raster; where bounding box shapefile is a polygon for raster layer extent (useful to "visualize" the complete raster).

Now, the issue is how "to put" a point in each raster cell corresponding to its perimeter. Next code, with methods of GDAL Python module, it does that work. Part implemented with PyQGIS was only for visualization purpose (because this question is for ArcGis).

from osgeo import gdal

dataset = gdal.Open('/home/zeito/pyqgis_data/aleatorio_mask.tif')

band = dataset.GetRasterBand(1)
noData = band.GetNoDataValue()

print noData

cols = band.XSize
rows = band.YSize

print rows, cols

transform = dataset.GetGeoTransform()

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

data = band.ReadAsArray(0, 0, cols, rows)

idxs = []

try:
    for i, row in enumerate(data):
        for j, element in enumerate(row):
            if data[i-1][j] == noData and element != noData or \
               data[i+1][j] == noData and element != noData or \
               data[i][j-1] == noData and element != noData or \
               element != noData and data[i][j+1] == noData:
               idxs.append([i, j])

except IndexError:
    print "Done!"

x0 = xOrigin + pixelWidth/2
y0 = yOrigin - pixelHeight/2

points = []

for item in idxs:
    x = x0 + item[1]*pixelWidth
    y = y0 - item[0]* pixelHeight
    points.append(QgsPoint(x,y))

wkt = dataset.GetProjection()

crs = QgsCoordinateReferenceSystem()
crs.createFromWkt(wkt)

epsg = crs.postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)   

After running the code at the Python Console of QGIS, I got:

Now, points list (associated to its raster values) in my code can be used for finding nearest raster cell value, for whatever vector point layer, considering: a) a point over nodata values b) a point out of raster extent.