Skip to content

Instantly share code, notes, and snippets.

@denis-bz

denis-bz/0-warp_raster_test.md

Last active September 13, 2021 11:02
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save denis-bz/0279a280a3221b1a96a7276290a5dd95 to your computer and use it in GitHub Desktop.
Save denis-bz/0279a280a3221b1a96a7276290a5dd95 to your computer and use it in GitHub Desktop.
Download ZIP
rioxarray reproject UTM -> 4326 widens NaNs at the sides 13 Sep 2021
Raw
0-warp_raster_test.md

rioxarray reproject UTM -> 4326 widens NaNs at the sides

Keywords: rioxarray, rasterio, reproject, geotiff, testcase

This plot shows that rioxarray reproject UTM -> 4326, middle -> right, widens NaNs at the sides:

13sep-warp_raster_test_riox

NaN widening also shows up on geotiffs from globalwindatlas.

I may well have done something stupid here -- new to gis. Reproject has lots of options in rioxarray / rasterio / ... don't know where this bottoms out. Do any of these have

  • a narrative description that one can follow with pencil and paper
  • a reference implementation for simple cases, using only numpy and scipy
  • a thorough test suite, not just eyeball plots ?

(Roundtrip tests are easy and I think useful, but not relevant here; I started with a 4326 raster just because it's easy and standalone.)

Code and .log is under my gists.

show_versions

rioxarray (0.5.0) deps:
  rasterio: 1.2.6
    xarray: 0.19.0
      GDAL: 3.3.0

Other python deps:
     scipy: 1.7.1
    pyproj: 3.1.0

System:
    python: 3.7.9

cheers
-- denis-bz-py@t-online.de 13 Sept 2021

Raw
rioutil.py
#!/usr/bin/env python3
""" rioutil.py: gen_xarray printxds ... rubble """
from itertools import product
import re
import numpy as np
from pyproj import CRS # https://pyproj4.github.io/pyproj/stable/examples.html
import rioxarray # https://corteva.github.io/rioxarray
from xarray import DataArray # here ~ namedtuple( y x values crs )
_lat = np.arange( 72, 60 - .01 )
_lon = np.arange( 20, 30.01 )
def gen_xarray( y=_lat, x=_lon, vals=None, crs=4326, ymul=0, verbose=1) -> DataArray:
y = np.asarray( y ) # dtype=float ?
x = np.asarray( x )
if vals is None:
vals = np.add.outer( ymul * y, x ) # ymul 100: 6020 6021 .. 6120 6121 ..
xa = DataArray( vals, coords=dict( y=y, x=x ))
if crs == 0:
xa.rio.set_crs( 4326 )
crs = xa.rio.estimate_utm_crs()
xa = xa.rio.write_crs( crs ) # set_crs ?
# https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html
xa = xa.rio.write_nodata ( np.NaN )
if verbose:
printxds( xa, "gen_xarray", verbose=verbose )
return xa
def printxds( xa: DataArray, txt="", verbose=2 ):
""" print a few summary lines of a geo xarray """
x = xa.x.values # may be nonuniform / jitter
y = xa.y.values
vals = xa.values
print( "\n{", txt, vals.shape )
print( "y: %.8g .. %.8g Δ %.8g " % (
y[0], y[-1], y[1] - y[0] ))
print( "x: %.8g .. %.8g Δ %.8g " % (
x[0], x[-1], x[1] - x[0] ))
print( "values:", vals.shape, vals.dtype )
print( " ", _quantiles( vals ))
if verbose:
print( vals ) # caller's printoptions
percentnan = np.isnan( vals ).sum() / vals.size * 100
if percentnan > 0:
print( "nan: %.1f %%" % percentnan )
print( "bounds w s e n: %g %g %g %g " % xa.rio.bounds() )
crs = xa.rio.crs
if verbose >= 2:
print( crs_wkt( crs ))
else:
print( "crs:", crs )
print( repr( xa.rio.transform() ))
# Affine(0.002500000000000124, 0.0, 3.349072682193763,
# 0.0, -0.002500000000000124, 55.91826874773593)
print( "}\n" )
def read_tif( geotiff: "filename", verbose=1 ) -> DataArray:
# geotiff = nu.findfile( geotiff, dirs= )
xa = rioxarray.open_rasterio( geotiff, parse_coordinates=True ).squeeze()
if verbose:
printxds( xa, geotiff )
return xa
# $rioxarray/_io.py xy shifted to pixel centers, -> Pixel is area
def write_tif( tifout, xa: DataArray ):
xa.rio.to_raster( tifout )
def crs_wkt( epsg ):
crs = CRS.from_user_input( epsg ) # int dict str
wkt = crs.to_wkt( pretty=True )
# rasterio.crs.CRS.from_epsg( epsg ) no pretty
wkt = re.sub( r"\n^ *[IOL]", lambda m: " " + m[0].strip(), wkt, flags=re.M )
return "{ %s \n%s \n}" % (epsg, wkt)
def utm_crs( zone: int ):
# 32 -> 32600 + 632 (32N == 32V ?)
crs = CRS.from_string( "+proj=utm +zone=%d +north" % zone )
return crs, crs.to_authority()[1]
def clip_yyxx( xa: DataArray, y0, y1, x0, x1: float ):
# box primitives intersect union ... ?
# x0, x1 = xa.x.values.min(), xa.x.values.max()
# y0, y1 = yds.y.values.min(), yds.y.values.max()
return xa.rio.clip_box(
miny=y0, maxy=y1, minx=x0, maxx=x1 )
def clip_pq( xa: DataArray, p, q: "Pt" ):
return xa.rio.clip_box(
miny=p[0], minx=p[1], maxy=q[0], maxx=q[1] )
def km( xa: DataArray ):
return xa.assign_coords( dict(
x = (xa.x - 500000) / 1000,
y = xa.y / 1000 ))
# set units km ?
def trim_nan( xa: DataArray, verbose=1 ):
""" trim all-NaN border rows and columns """
A = xa.values.squeeze()
assert len( A.shape ) == 2, A.shape
Anan = np.isnan( A )
nanrows = np.all( Anan, axis=1 ) # fast first and last ?
jrow0, jrow1 = (~ nanrows) .nonzero()[0] [[0, -1]]
nancols = np.all( Anan, axis=0 )
jcol0, jcol1 = (~ nancols) .nonzero()[0] [[0, -1]]
xtrim = xa[ jrow0 : jrow1+1, jcol0 : jcol1 + 1 ] # with attrs
if verbose:
print( "\ntrim_nan: %s -> %s y %s -> %s x %s -> %s \n" % (
A.shape, xtrim.shape,
nu.lohif( xa.y ), nu.lohif( xtrim.y ),
nu.lohif( xa.x ), nu.lohif( xtrim.x )))
return xtrim
#...............................................................................
# etc --
def linstep( a, b, step=1 ) -> "linspace including b":
n = round( abs( (b - a) / step )) + 1
return np.linspace( a, b, int( n ))
def allpairs( *arrays ):
return np.array( list( product( *arrays )))
# allpairs( [0,1,2], [3,4] ) # [[0 3] [0 4] [1 3] [1 4] [2 3] [2 4]]]
# return np.stack( np.meshgrid(*arrays), axis=-1 ).reshape(-1, len(arrays))
# allpairs( [0,1,2], [3,4] ) # [[0 3] [1 3] [2 3] [0 4] [1 4] [2 4]]
# https://stackoverflow.com/questions/11144513/cartesian-product-of-x-and-y-array-points-into-single-array-of-2d-points
_q = np.arange( 0, 100.01, 10 )
def _quantilesf( X, q=_q, ints=0 ):
""" np.percentile opt astype int """
p = np.nanpercentile( X, q=q )
if ints != 0:
p = (p * ints) .round().astype(int)
return p
def _quantiles( X, q=_q, ints=0 ):
ints = np.nanmax( X ) < 1e4
return "quantiles: %s" % _quantilesf( X, q=q, ints=ints )
Raw
warp_raster_test_riox-13sep.log
# from: warp_raster_test_riox.py
# 13 Sep 2021 10:35z py/geo/rio Denis-iMac 10.10.5
▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄
rioxarray (0.5.0) deps:
rasterio: 1.2.6
xarray: 0.19.0
GDAL: 3.3.0
Other python deps:
scipy: 1.7.1
pyproj: 3.1.0
System:
python: 3.7.9 (v3.7.9:13c94747c7, Aug 15 2020, 01:31:08) [Clang 6.0 (clang-600.0.57)]
executable: /Library/Frameworks/Python.framework/Versions/3.7/bin/python
machine: Darwin-14.5.0-x86_64-i386-64bit
warp_raster_test_riox.py lat (72, 60) lon (20, 30) step 1 zoom None Resampling.bilinear
lat: [72 71 70 69 68 67 66 65 64 63 62 61 60]
lon: [20 21 22 23 24 25 26 27 28 29 30]
{ gen_xarray (13, 11)
y: 72 .. 60 Δ -1
x: 20 .. 30 Δ 1
values: (13, 11) float64
quantiles: [20 21 22 23 24 25 26 27 28 29 30]
[[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]
[20 21 22 23 24 25 26 27 28 29 30]]
bounds w s e n: 19.5 59.5 30.5 72.5
crs: EPSG:4326
Affine(1.0, 0.0, 19.5,
0.0, -1.0, 72.5)
}
{ reproject EPSG:32635 Resampling.bilinear None (16, 7)
y: 8015571.9 .. 6670511.3 Δ -89670.712
x: 120796.55 .. 658820.82 Δ 89670.712
values: (16, 7) float64
quantiles: [20 21 22 23 24 25 26 27 28 29 30]
[[nan nan 21.1 23.8 26.3 29 nan]
[nan nan 21.3 23.9 26.4 29 nan]
[nan nan 21.7 24 26.4 28.9 nan]
[nan 20.1 21.9 24.1 26.4 28.8 nan]
[nan 20.2 22 24.2 26.5 28.8 nan]
[nan 20.3 22.1 24.3 26.5 28.7 nan]
[nan 20.4 22.3 24.4 26.5 28.7 nan]
[nan 20.5 22.5 24.5 26.6 28.6 nan]
[nan 20.7 22.7 24.6 26.6 28.5 29.9]
[nan 20.9 22.8 24.7 26.6 28.4 29.9]
[nan 21.1 22.9 24.8 26.6 28.4 29.8]
[nan 21.2 23 24.9 26.7 28.3 29.8]
[20.2 21.3 23.1 24.9 26.7 28.3 29.8]
[20.2 21.6 23.2 24.9 26.7 28.3 29.7]
[20.3 21.7 23.3 25 26.7 28.2 29.7]
[20.4 21.8 23.4 25 26.7 28.2 29.7]]
nan: 20.5 %
bounds w s e n: 75961.2 6.62568e+06 703656 8.06041e+06
{ EPSG:32635
PROJCRS["WGS 84 / UTM zone 35N",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]],
CONVERSION["UTM zone 35N",
METHOD["Transverse Mercator", ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",27,
ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1], ID["EPSG",8805]],
PARAMETER["False easting",500000, LENGTHUNIT["metre",1], ID["EPSG",8806]],
PARAMETER["False northing",0, LENGTHUNIT["metre",1], ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1]],
AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1]], ID["EPSG",32635]]
}
Affine(89670.71194328342, 0.0, 75961.18991876603,
0.0, -89670.71194328342, 8060407.292590056)
}
{ reproject back to EPSG:4326 (11, 16)
y: 72.054904 .. 60.348155 Δ -1.1706749
x: 15.046899 .. 32.607022 Δ 1.1706749
values: (11, 16) float64
quantiles: [20 21 22 23 24 24 26 27 28 29 30]
[[nan nan nan nan nan 21.2 22.1 23.3 24.4 25.6 26.7 27.9 29 29 nan nan]
[nan nan nan nan nan 21.5 22.1 23.2 24.4 25.6 26.7 27.9 28.9 nan nan nan]
[nan nan nan 20.1 20.3 21.2 22.1 23.2 24.4 25.6 26.7 27.9 28.8 nan nan nan]
[nan nan nan nan 20.2 21.1 22.1 23.2 24.4 25.6 26.8 28 28.8 nan nan nan]
[nan nan nan nan 20.4 21 22.1 23.2 24.4 25.6 26.8 28 28.7 nan nan nan]
[nan nan nan nan 20.6 20.9 22.1 23.3 24.4 25.6 26.8 28 28.8 nan nan nan]
[nan nan nan nan nan 20.9 22.1 23.3 24.5 25.6 26.8 27.9 28.9 29.8 nan nan]
[nan nan nan nan nan 21.1 22.1 23.3 24.5 25.6 26.8 27.9 28.9 29.8 nan nan]
[nan nan nan nan 20.3 21.1 22.1 23.2 24.4 25.6 26.8 27.9 28.9 29.8 nan nan]
[nan nan nan nan 20.2 21.1 22.1 23.2 24.4 25.6 26.8 27.9 28.9 29.7 nan nan]
[nan nan nan nan 20.3 21 22.1 23.2 24.4 25.6 26.8 27.9 29 29.7 nan nan]]
nan: 42.0 %
bounds w s e n: 14.4616 59.7628 33.1924 72.6402
{ EPSG:4326
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north, ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east, ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]]
}
Affine(1.1706748886863299, 0.0, 14.46156130403454,
0.0, -1.1706748886863299, 72.64024118457897)
}
Raw
warp_raster_test_riox.py
#!/usr/bin/env python3
r""" warp_rester_test_riox roundtrip: %s $\, \to \,$ %s $\, \to \,$ %s """
import numpy as np
import matplotlib.pyplot as pl
import seaborn as sns
import rioxarray
# https://corteva.github.io/rioxarray/stable/_modules/rioxarray/raster_array.html#RasterArray.reproject
# -> rasterio.warp.reproject
from rasterio.enums import Resampling
from xarray import DataArray # here just a type
from rioutil import gen_xarray, linstep, printxds
#...............................................................................
def reproject2_riox( xa: DataArray, tocrs,
resampling=Resampling.bilinear, shape=None ):
""" -> proj, proj2
proj = reproject( xa )
proj2 = reproject( proj ) -- roundtrip test
"""
crs = xa.rio.crs
if tocrs == 0:
tocrs = xa.rio.estimate_utm_crs()
proj = xa.rio.reproject( tocrs, resampling=resampling, shape=shape )
# transform= ?
printxds( proj, "reproject %s %s %s " % (
tocrs, resampling, shape ))
proj2 = proj.rio.reproject( crs, resampling=resampling )
printxds( proj2, "reproject back to %s" % crs )
return proj, proj2
#...............................................................................
if __name__ == "__main__":
import sys
np.set_printoptions( edgeitems=10, threshold=10, linewidth=130,
formatter = dict( float = lambda _: "%.3g" % _ ))
thispy = sys.argv[0]
print( 80 * "▄" )
# print( "▄ bp" ); breakpoint()
rioxarray.show_versions()
#...........................................................................
lat = 72, 60 # Finland ~ 58.8 .. 70, 19 .. 31.6
lon = 20, 30
step = 1
zoom = None
tocrs = 0 # 32635
resampling = Resampling.bilinear # nearest bilinear cubic cubic_spline no diff
plot = 0
save = 0
vmin = vmax = None
# to change these params, run this.py a=1 b=None 'c = expr' ... in sh or ipython --
for arg in sys.argv[1:]:
exec( arg )
print( "%s lat %s lon %s step %g zoom %s %s " % (
thispy, lat, lon, step, zoom, resampling ))
lat = linstep( *lat, step ) # linspace [lo .. hi]
lon = linstep( *lon, step )
print( "lat:", lat )
print( "lon:", lon )
#...........................................................................
xa = gen_xarray( lat, lon, crs=4326 ) # lat x lon rioxarray
if zoom is not None:
zoom = (zoom * xa.shape[0] - 1, zoom * xa.shape[1] - 1)
proj, proj2 = reproject2_riox( xa, tocrs=tocrs, resampling=resampling, shape=zoom )
crs = xa.rio.crs
tocrs = proj.rio.crs
if plot:
fig, (ax0, ax1, ax2) = pl.subplots( ncols=3, figsize=[11.7, 8.3] )
title = __doc__ % (crs, tocrs, crs)
fig.suptitle( title.strip(), multialignment="left" )
kw = dict( cmap="plasma", vmin=vmin, vmax=vmax )
xa.plot( ax=ax0, **kw ) # drop colorbar ?
ax0.set_title( "%s %s" % (crs, xa.shape ))
proj.plot( ax=ax1, **kw )
ax1.set_title( r"to %s %s, north narrower" % (proj.rio.crs, proj.shape ))
proj2.plot( ax=ax2, **kw )
ax2.set_title( r"to %s %s" % (crs, proj2.shape ))
if plot >= 2:
import zutil as nu
nu.savefig( __file__ )
if save:
print( "\n> tmp.tif tmp-proj.tif tmp-proj2.tif" )
xa.rio.to_raster( "tmp.tif" )
proj.rio.to_raster( "tmp-proj.tif" )
proj2.rio.to_raster( "tmp-proj2.tif" )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment