-
Notifications
You must be signed in to change notification settings - Fork 551
Description
Expected behavior and actual behavior.
I've wrapped rasterio.wrap.reproject to work on two-dimensional xarray DataArrays, and I've written some unit tests. My CI started breaking yesterday, and I can trace it back to the upgrade from rasterio 1.0.28 to 1.1.0.
At first I thought the break was caused by a change in rasterio.warp.reproject, but it looks like it's the changes in constructing the CRS instead. It looks like these changes are changes for the better (fortunately), but I don't think they show up in the changelog. So this issue is primarily for my own understanding, but possibly highlights an unaware change.
I've installed two conda environments, one with 1.0.28 and one with 1.10.0, which result in the examples below.
Steps to reproduce the problem.
The first change is this in the from_string method:
import rasterio
rasterio.crs.CRS.from_string("+init=epsg:28992")On 1.0.28 this returns:
CRS.from_epsg(28992)On 1.10.0 this returns:
CRS.from_wkt('PROJCS["Amersfoort / RD New",GEOGCS["Amersfoort",DATUM["Amersfoort",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6289"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4289"]],PROJECTION["Oblique_Stereographic"],PARAMETER["latitude_of_origin",52.1561605555556],PARAMETER["central_meridian",5.38763888888889],PARAMETER["scale_factor",0.9999079],PARAMETER["false_easting",155000],PARAMETER["false_northing",463000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')There's also the deprecation warning:
Warning 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.Anyway, so far so good.
However, when using the crs object to reproject, the results are different.
import affine
import numpy as np
import rasterio
import rasterio.warp
nrow = 5
ncol = 8
source = np.arange(nrow * ncol).reshape(nrow, ncol)
src_crs = rasterio.crs.CRS.from_string("+init=epsg:28992")
src_transform = affine.Affine(1.0, 0.0, 20.0, 0.0, -1.0, 30.0)
dst_crs = {"init": "EPSG:32631"}
src_height = nrow
src_width = ncol
bounds = rasterio.transform.array_bounds(src_height, src_width, src_transform)
dst_transform, dst_width, dst_height = rasterio.warp.calculate_default_transform(
src_crs, dst_crs, src_width, src_height, *bounds
)
destination = np.empty((nrow * 2, ncol * 2))
rasterio.warp.reproject(
source=source,
destination=destination,
src_transform=src_transform,
dst_transform=dst_transform,
src_crs=src_crs,
dst_crs=dst_crs,
resampling=rasterio.enums.Resampling.nearest,
)
print(destination)In 1.0.28:
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
In 1.10.0:
Warning 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
[[ 0. 1. 2. 3. 4. 5. 6. 7. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 9. 10. 11. 12. 13. 14. 15. 0. 0. 0. 0. 0. 0. 0. 0.]
[16. 17. 18. 19. 20. 21. 22. 23. 0. 0. 0. 0. 0. 0. 0. 0.]
[24. 25. 26. 27. 28. 29. 30. 31. 0. 0. 0. 0. 0. 0. 0. 0.]
[32. 33. 34. 35. 36. 37. 38. 39. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
Basically, it appears to be failing to reproject in 1.0.28 when creating the crs from the string, and it works in 1.10.0.
In 1.0.28, the REPL representation is crs.from_epsg(28992). So, what if I force both to from_espg?
In the scrip above, redefining src_crs:
#src_crs = rasterio.crs.CRS.from_string("+init=epsg:28992")
src_crs = rasterio.crs.CRS(init="epsg:28992")REPL says CRS.from_epsg(28992) both times, and I get identical results, with non-zero values both times.
However, in both 1.0.28 and 1.10.0, this evaluates as True:
rasterio.crs.CRS(init="epsg:28992") == rasterio.crs.CRS.from_string("+init=epsg:28992")So I'm really not sure what's going on. It looks like a bug has been fixed (since the array wasn't projected properly in 1.0.28 -- one on which my CI was (incorrectly!) relying.)
Operating system
OS: Windows 10 64-bit Enterprise, Version 1803
Rasterio version and provenance
rasterio versions 1.0.28 and 1.10.0 installed using conda install -c conda-forge