python example for the basics

[mdsumner]

https://gist.github.com/vincentsarago/b00f50f8b66ab5ccbd4449f6a1bd8b71#file-rasterize_point-py-L79

from this tweet (they missed the IDW part, but still very useful code for me) https://twitter.com/_VincentS_/status/1597971755939016710?s=20&t=ZpcAHSd0zMVO8K0iAZjq0Q

"""Rasterize Points

requirements:
- rasterio
- sparse

"""

import typing
import math

import warnings

import numpy

import rasterio
from rasterio.crs import CRS
from rasterio.dtypes import _gdal_typename
from rasterio.enums import Resampling
from rasterio.shutil import copy
from rasterio.transform import Affine, from_bounds
from rasterio.warp import transform as transform_points
from sparse import COO

in_crs = CRS.from_epsg(4326)
out_crs = CRS.from_epsg(3857)

# Point values in form of value, lon, lat
point_values: typing.List[typing.Tuple[float, float, float]]

# unpack points
values, lon, lat = zip(*point_values)

# Translate the lat/lon coordinates to `out_crs`
x, y = transform_points(in_crs, out_crs, lon, lat)

nodata_value: typing.Any = 0

# Either define the output resolution (in out_crs projection) or the output shape
out_shape: typing.Optional[typing.Tuple[int, int]]
resolution: typing.Optional[float]

xmin = numpy.amin(x)
xmax = numpy.amax(x)
ymin = numpy.amin(y)
ymax = numpy.amax(y)
bounds = [xmin, ymin, xmax, ymax]

if out_shape is not None:
    height, width = out_shape
    resolution = max(
        (bounds[2] - bounds[0]) / width,
        (bounds[3] - bounds[1]) / height,
    )

elif resolution is not None:
    width = round((bounds[2] - bounds[0]) / resolution) + 1
    height = round((bounds[3] - bounds[1]) / resolution) + 1

else:
    raise Exception("missing output resolution or shape")

# Output Affine Transform
transform = Affine.translation(bounds[0] - resolution / 2, bounds[1] - resolution / 2) * Affine.scale(resolution, resolution)

# Get Row/Col indexes for each point
rows, cols = rasterio.transform.rowcol(transform, x, y, op=math.floor)  # return rows, cols

# Rasterize Dense Matrix using Sparse COO
data = COO((rows, cols), values, fill_value=nodata_value).todense()

# Make sure width/height match array shapes
assert (height, width) == data.shape

# Create the Output Raster
info = {
    "DATAPOINTER": data.__array_interface__["data"][0],
    "PIXELS": width,
    "LINES": height,
    "BANDS": 1,
    "DATATYPE": _gdal_typename(data.dtype.name),
}

strides = data.__array_interface__.get("strides", None)
if strides is not None:
    if len(strides) == 2:
        lineoffset, pixeloffset = strides
        info.update(LINEOFFSET=lineoffset, PIXELOFFSET=pixeloffset)
    else:
        bandoffset, lineoffset, pixeloffset = strides
        info.update(
            BANDOFFSET=bandoffset,
            LINEOFFSET=lineoffset,
            PIXELOFFSET=pixeloffset,
        )

dataset_options = ",".join(f"{name}={val}" for name, val in info.items())
datasetname = f"MEM:::{dataset_options}"
with warnings.catch_warnings():
    warnings.simplefilter("ignore", rasterio.errors.NotGeoreferencedWarning)
    with rasterio.open(datasetname, "r+") as src:
        src.crs = out_crs
        src.transform = transform
        src.nodata = nodata_value

        dst_profile = {
            "driver": "COG",
            "blocksize": 256,
            "compress": "DEFLATE",
            "sparse_ok": "TRUE",
            "nodata": nodata_value,
        }
        fout = f"cog.tif"
        copy(src, fout, **dst_profile)

[mdsumner]

that code does this

n <- 1000
pts0 <- cbind(runif(n, -180, 180), runif(n, -80, 80))
in_crs <-  "OGC:CRS84"
out_crs <- "EPSG:3857"

pts <- reproj::reproj_xy(pts0, out_crs, source = in_crs)

## grid specification (bounds, shape)
## we could be tidier about the extent, but I'm not sure how the py was intended to user-input this
ex <- c(range(pts[,1]), range(pts[,2]))
dm <- c(200, 200)  ## again, could be aspect ratioed of the data but just go for it

## the rasterio transform point to grid part
cell <- vaster::cell_from_xy(dm, ex, pts)

## we're not sparse now
cnt <- matrix(tabulate(cell, nbins = prod(dm)), dm[2L], byrow = TRUE)

## now we need something to write a tif
vapour::vapour_create("file.tif", extent = ex, dimension = dm, projection = out_crs, n_bands = 1, overwrite = TRUE)
## I think the current orientation is the write fun might be a mistake ... or at least needs doc
vapour::vapour_write_raster_block("file.tif", t(cnt), offset = c(0, 0), dimension = dm, band = 1, overwrite = TRUE)

## now, get the info and plot
info <- vapour::vapour_raster_info("file.tif")
dat <- whatarelief::elevation(source = "file.tif", extent = info$extent, dimension = info$dimension, projection = info$projection)
ximage::ximage(dat, extent = info$extent, asp = 1)
points(pts)