denis-bz/0-warp.md

## 0-warp.md

      
    Raw
  

              0-warp.md
            
          
    Warping maps or images

Keywords: image raster resize rescale reproject Manhattan-grid regular-grid
We have a picture (image, raster) A, an array of colored pixels,
and want to warp it to a new picture O of a different size or shape.
Examples:
Map projection, flat maps <-> globe;
Image scaling. 

There are hundreds of cases of these, and zillions of programs
with a huge range of capabilities, styles, and learning curves.
The purpose of this note is to describe a bare-bones warp.py for students,
with pictures A and O on Manhattan grids -- simpler than general grids.
An example of a round-trip test, 4326 -> utm31 -> 4326, with warp.py:

Maps and grids, snaptogrid

Think of a map on two layers:

a flat paper map with continuous coordinates,
say degrees or metres on a map of the earth's surface,
or millimetres on a computer screen.
a Manhattan grid: a transparent plastic sheet overlaying that,
with streets running east-west at y[0] y[1] ...
and avenues north-south at x[0] x[1] ... .
x and y may be uniformly spaced, np.linspace, or may not be.

Given a bunch of points on a map / in Manhatten, there's an important function
snaptogrid: each point -> the nearest streetcorner on the grid, indices `i j`.

O grid -> A grid

We have pixels on the whole A grid, and warp transforms both ways, AtoO and OtoA.
(These may be nonlinear, but are very nearly inverses --
AtoO( OtoA( points )) ~= the same points.) 

How do we warp the pixels on the A grid to an O grid
of different size or shape ?
One way is this pipeline, in warp.py:

transform O gridpoints -> a cloud in A
| trim points ouside A.bbox
| snap to A gridpoints -- maybe dups, maybe missing
| Apix[I,J] at A gridpoints
| Opix at the corresponding O gridpoints


Types

Type names may help to keep track of what's what (or may be mttiw).
Most things in this implementation are numpy ndarrays, 1d or n × 2,
of various kinds:
Pts = "nx2 floats"  # a common idiom: y, x = yx.T
Gridpts = Pts   on a Manhattan grid (street corners),
                Y in y[0] y[1] .. y[ny-1], X in x[0] x[1] .. x[nx-1]
Cloud = Pts     anywhere -- on Gridpts, between them, or outside the bounding box
Ij = "nx2 int indices"   I, J = ij.T
Gridarray = "2d numpy array"  values A[I,J] at Gridpts

snaptogrid: "Cloud -> Ij"
Transform = "Cloud -> Cloud"  may be nonlinear, may go outside a bounding box

Notes

In pseudocode on paper, it's useful to think of arrays indexed by continuous y x.
At the end of the day, real code will have to map floats y x on the continuous layer
<-> int indices i j on the grid; these maps are simple but really gum up code.
There are many ways to implement snaptogrid.
warp.py uses Near_rgrid, just a page of code;
scipy RegularGridInterpolator;
also does bilinear interpolation;
scipy griddata ==
scipy KDTree
does points -> k nearest neighbors in an arbitrary point cloud.
See also: 

rioxarray reproject UTM -> 4326 widens NaNs at the sides
-- bug, or feature ?
Comments welcome

will trade the complete code, about 500 lines, for test cases.
cheers 

-- denis-bz-py@t-online.de  1 Nov 2021

  
## 1-warp.py
#!/usr/bin/env python3
""" warp( A grid -> O grid, transformOtoA, Apix ), see 0-warp.md """
# keywords: image raster resize rescale reproject Manhattan-grid regular-grid python

import numpy as np

from rgrid import Rgrid  # Manhattan aka regular: .y .x .shape .bbox .snaptogrid
from bbox import inbbox
from util import allpairs, allpairsij


def warp(
    A: "Manhattan grid in: pairs y_i x_j",
    O: "Manhattan grid out",
    transformOtoA: "Gridpoints in O -> Cloud in A",  # e.g. pyproj.Transformer.from_crs
    Apix: "2d ndarray of pixels at A gridpoints",
    verbose=1
    ) -> "2d ndarray of pixels at O gridpoints":

    """ pipeline:
    transform O gridpoints -> a cloud in A
    | trim points ouside A.bbox
    | snap to A gridpoints, streetcorners in Manhattan -- maybe dups, maybe missing
    | Apix[I,J] at A gridpoints
    | Opix at the corresponding O gridpoints
    """
    assert A.shape == Apix.shape, [A.shape, Apix.shape]

    oyx = allpairs( O.y, O.x )
    oy, ox = oyx.T
    oij = allpairsij( *O.shape )

    #...........................................................................
    ayx = transformOtoA( oy, ox )  # may be nonlinear

    ain = inbbox( ayx, A.bbox )
    pcin = len(ain) / len(ayx) * 100
    if pcin < 100:
        ayx = ayx[ain]
        oij = oij[ain]

    aI, aJ = A.snaptogrid( ayx )  # nearest streetcorners in Manhattan
    apix = Apix[ aI, aJ ]  # flat

    opix = np.full( O.shape, np.NaN )
    oI, oJ = oij.T
    opix[ oI, oJ ] = apix
    if verbose:
        print( "\nwarp: A %s  O %s  %s  %.0f %% inside A.bbox" % (
                A.shape, O.shape, transformOtoA.__name__, pcin ))

    return opix


__version__ = "2021-11-02-Nov"  # denis-bz-py t-online.de


## near_rgrid.py
#!/usr/bin/env python3
# keywords: nearest-neighbor Manhattan-grid regular-grid python numpy Voronoi

import numpy as np

#...............................................................................
class Near_rgrid( object ):
    """ nearest neighbors on a Manhattan aka regular grid
    1d:
    near = Near_rgrid( x: sorted 1d array )
    nearix = near.query( q: 1d ) -> indices of the points x_i nearest each q_i
        x[nearix[0]] is the nearest to q[0]
        x[nearix[1]] is the nearest to q[1] ...
        nearpoints = x[nearix] is near q
    If A is an array of pixels / colors at x[0] x[1] ...,
        A[nearix] are the colors near q[0] q[1] ...

    2d: on a Manhattan aka regular grid,
        streets running east-west at y_i, avenues north-south at x_j,
        we want to find the streetcorners nearest some query points:
    near = Near_rgrid( y, x: sorted 1d arrays, e.g. latitide longitude )
    I, J = near.query( q: nq × 2 array, columns qy qx )
    -> nq × 2 indices of the nearest streetcorners
        gridpoints = np.column_stack(( y[I], x[J] ))  # e.g. street corners
        diff = gridpoints - querypoints
        distances = norm( diff, axis=1, ord= )
    If A is an array of colors at gridpoints,
        A[I,J] are the colors at grid points nearest the query points q[0] q[1] ...

    near.snap is an alias for near.query, [[y x] ...] -> indices [[i j] ...]
    near.snaptoyx -> grid points [[y_i x_j] ...]
    Queries outside a grid snap to the nearest on-grid, e.g. Bronx -> Manhattan.

    See Howitworks below, and the plot Voronoi-random-regular-grid.
    """

    def __init__( self, *axes: "1d arrays" ):
        axarrays = []
        for ax in axes:
            axarray = np.asarray( ax ).squeeze()
            assert axarray.ndim == 1, "each axis should be 1d, not %s " % (
                    str( axarray.shape ))
            axarrays += [axarray]
        self.midpoints = [_midpoints( ax ) for ax in axarrays]
        self.axes = axarrays
        self.ndim = len(axes)

    def query( self, queries: "nq × dim points" ) -> "nq × dim grid indices":
        """ -> the indices of the nearest points in the grid """
        queries = np.asarray( queries ).squeeze()  # or list x y z ?
        if self.ndim == 1:
            assert queries.ndim <= 1, queries.shape
            return np.searchsorted( self.midpoints[0], queries )  # scalar, 0d ?
        queries = np.atleast_2d( queries )
        assert queries.shape[1] == self.ndim, [
                queries.shape, self.ndim]
        return [np.searchsorted( mid, q )  # parallel: k axes, k processors
                for mid, q in zip( self.midpoints, queries.T )]

    snaptogrid = ijnear = query

    def yxnear( self, queries: "nq × dim points" ) -> "nq × dim grid points":
        """ -> the nearest points in the grid, 2d [[y_j x_i] ...] """
        ix = self.query( queries )
        if self.ndim == 1:
            return self.axes[0][ix]
        else:
            axix = [ax[j] for ax, j in zip( self.axes, ix )]
            return np.array( axix ).squeeze()


def _midpoints( points: "array-like 1d, *must be sorted*" ) -> "1d":
    points = np.asarray( points ).squeeze()
    assert points.ndim == 1, points.shape
    diffs = np.diff( points )
    mindiff = np.nanmin( diffs )
    assert mindiff > 0, \
        "Near_rgrid: the input array must be sorted, not %g .. %g  min diff %.2g " % (
        points[0], points[-1], mindiff )
    return (points[:-1] + points[1:]) / 2  # floats

#...............................................................................
Howitworks = \
"""
How Near_rgrid works in 1d:
Consider the midpoints halfway between fenceposts | | |
The interval [left midpoint .. | .. right midpoint] is what's nearest each post --

    |   |       |                     |   points
    | . |   .   |          .          |   midpoints
      ^^^^^^               .            nearest points[1]
            ^^^^^^^^^^^^^^^             nearest points[2]  etc.

2d:
    I, J = Near_rgrid( y, x ).query( q )
    I = nearest in `y`, J = nearest in `x`
    A[I,J]: pixels in array A, near the query points.

Notes
-----
If a query point is exactly halfway between two data points,
e.g. on a grid of ints, the lines (x + 1/2) U (y + 1/2),
which "nearest" you get is implementation-dependent, unpredictable.

Nearest in any norm ? See the plot Voronoi-random-regular-grid.

"""

Murky = \
""" NaNs in points, in queries ?
"""

__version__ = "2021-10-30 oct  denis-bz-py"
	#!/usr/bin/env python3
	""" warp( A grid -> O grid, transformOtoA, Apix ), see 0-warp.md """
	# keywords: image raster resize rescale reproject Manhattan-grid regular-grid python

	import numpy as np

	from rgrid import Rgrid # Manhattan aka regular: .y .x .shape .bbox .snaptogrid
	from bbox import inbbox
	from util import allpairs, allpairsij


	def warp(
	A: "Manhattan grid in: pairs y_i x_j",
	O: "Manhattan grid out",
	transformOtoA: "Gridpoints in O -> Cloud in A", # e.g. pyproj.Transformer.from_crs
	Apix: "2d ndarray of pixels at A gridpoints",
	verbose=1
	) -> "2d ndarray of pixels at O gridpoints":

	""" pipeline:
	transform O gridpoints -> a cloud in A
	\| trim points ouside A.bbox
	\| snap to A gridpoints, streetcorners in Manhattan -- maybe dups, maybe missing
	\| Apix[I,J] at A gridpoints
	\| Opix at the corresponding O gridpoints
	"""
	assert A.shape == Apix.shape, [A.shape, Apix.shape]

	oyx = allpairs( O.y, O.x )
	oy, ox = oyx.T
	oij = allpairsij( *O.shape )

	#...........................................................................
	ayx = transformOtoA( oy, ox ) # may be nonlinear

	ain = inbbox( ayx, A.bbox )
	pcin = len(ain) / len(ayx) * 100
	if pcin < 100:
	ayx = ayx[ain]
	oij = oij[ain]

	aI, aJ = A.snaptogrid( ayx ) # nearest streetcorners in Manhattan
	apix = Apix[ aI, aJ ] # flat

	opix = np.full( O.shape, np.NaN )
	oI, oJ = oij.T
	opix[ oI, oJ ] = apix
	if verbose:
	print( "\nwarp: A %s O %s %s %.0f %% inside A.bbox" % (
	A.shape, O.shape, transformOtoA.__name__, pcin ))

	return opix


	__version__ = "2021-11-02-Nov" # denis-bz-py t-online.de
	#!/usr/bin/env python3
	# keywords: nearest-neighbor Manhattan-grid regular-grid python numpy Voronoi

	import numpy as np

	#...............................................................................
	class Near_rgrid( object ):
	""" nearest neighbors on a Manhattan aka regular grid
	1d:
	near = Near_rgrid( x: sorted 1d array )
	nearix = near.query( q: 1d ) -> indices of the points x_i nearest each q_i
	x[nearix[0]] is the nearest to q[0]
	x[nearix[1]] is the nearest to q[1] ...
	nearpoints = x[nearix] is near q
	If A is an array of pixels / colors at x[0] x[1] ...,
	A[nearix] are the colors near q[0] q[1] ...

	2d: on a Manhattan aka regular grid,
	streets running east-west at y_i, avenues north-south at x_j,
	we want to find the streetcorners nearest some query points:
	near = Near_rgrid( y, x: sorted 1d arrays, e.g. latitide longitude )
	I, J = near.query( q: nq × 2 array, columns qy qx )
	-> nq × 2 indices of the nearest streetcorners
	gridpoints = np.column_stack(( y[I], x[J] )) # e.g. street corners
	diff = gridpoints - querypoints
	distances = norm( diff, axis=1, ord= )
	If A is an array of colors at gridpoints,
	A[I,J] are the colors at grid points nearest the query points q[0] q[1] ...

	near.snap is an alias for near.query, [[y x] ...] -> indices [[i j] ...]
	near.snaptoyx -> grid points [[y_i x_j] ...]
	Queries outside a grid snap to the nearest on-grid, e.g. Bronx -> Manhattan.

	See Howitworks below, and the plot Voronoi-random-regular-grid.
	"""

	def __init__( self, *axes: "1d arrays" ):
	axarrays = []
	for ax in axes:
	axarray = np.asarray( ax ).squeeze()
	assert axarray.ndim == 1, "each axis should be 1d, not %s " % (
	str( axarray.shape ))
	axarrays += [axarray]
	self.midpoints = [_midpoints( ax ) for ax in axarrays]
	self.axes = axarrays
	self.ndim = len(axes)

	def query( self, queries: "nq × dim points" ) -> "nq × dim grid indices":
	""" -> the indices of the nearest points in the grid """
	queries = np.asarray( queries ).squeeze() # or list x y z ?
	if self.ndim == 1:
	assert queries.ndim <= 1, queries.shape
	return np.searchsorted( self.midpoints[0], queries ) # scalar, 0d ?
	queries = np.atleast_2d( queries )
	assert queries.shape[1] == self.ndim, [
	queries.shape, self.ndim]
	return [np.searchsorted( mid, q ) # parallel: k axes, k processors
	for mid, q in zip( self.midpoints, queries.T )]

	snaptogrid = ijnear = query

	def yxnear( self, queries: "nq × dim points" ) -> "nq × dim grid points":
	""" -> the nearest points in the grid, 2d [[y_j x_i] ...] """
	ix = self.query( queries )
	if self.ndim == 1:
	return self.axes[0][ix]
	else:
	axix = [ax[j] for ax, j in zip( self.axes, ix )]
	return np.array( axix ).squeeze()


	def _midpoints( points: "array-like 1d, must be sorted" ) -> "1d":
	points = np.asarray( points ).squeeze()
	assert points.ndim == 1, points.shape
	diffs = np.diff( points )
	mindiff = np.nanmin( diffs )
	assert mindiff > 0, \
	"Near_rgrid: the input array must be sorted, not %g .. %g min diff %.2g " % (
	points[0], points[-1], mindiff )
	return (points[:-1] + points[1:]) / 2 # floats

	#...............................................................................
	Howitworks = \
	"""
	How Near_rgrid works in 1d:
	Consider the midpoints halfway between fenceposts \| \| \|
	The interval [left midpoint .. \| .. right midpoint] is what's nearest each post --

	\| \| \| \| points
	\| . \| . \| . \| midpoints
	^^^^^^ . nearest points[1]
	^^^^^^^^^^^^^^^ nearest points[2] etc.

	2d:
	I, J = Near_rgrid( y, x ).query( q )
	I = nearest in `y`, J = nearest in `x`
	A[I,J]: pixels in array A, near the query points.

	Notes
	-----
	If a query point is exactly halfway between two data points,
	e.g. on a grid of ints, the lines (x + 1/2) U (y + 1/2),
	which "nearest" you get is implementation-dependent, unpredictable.

	Nearest in any norm ? See the plot Voronoi-random-regular-grid.

	"""

	Murky = \
	""" NaNs in points, in queries ?
	"""

	__version__ = "2021-10-30 oct denis-bz-py"