corteva / rioxarray

geospatial xarray extension powered by rasterio
https://corteva.github.io/rioxarray
Other
530 stars 85 forks source link

QST: Reprojection not working properly #420

Closed g2giovanni closed 3 years ago

g2giovanni commented 3 years ago

Problem description

Hi, I'm not able to reproject properly a raster to EPSG:4326. As a result, I obtain a raster with (width, height) = (1, 1).

Following a snippet of code to reproduce the problem. You can download the raster for the test here: https://drive.google.com/file/d/1lcfKr3m17TnE1SOZOKCkplK8w0I55fZZ/view?usp=sharing

import rioxarray as rxr
import rasterio

INPUT_RASTER = 'aaaaaaa_64.tif'
OUPUT_RASTER = "aaaaaaa_64_4326.tif"

da = rxr.open_rasterio(INPUT_RASTER)

res = tuple(abs(res_val) for res_val in da[0].rio.resolution())
dest_crs=rasterio.crs.CRS.from_epsg(4326)

da.rio.reproject(dst_crs=dest_crs, resolution=res).rio.to_raster(OUPUT_RASTER)

Expected Output

Correct reprojected raster into EPSG:4326

Environment Information

snowman2 commented 3 years ago

The input raster has the projection EPSG:32632 which is in a projected coordinate system with units of meters.

import rioxarray

rds = rioxarray.open_rasterio("aaaaaaa_64.tif")
rds.rio.crs  # CRS.from_epsg(32632)
rds.rio.resolution()  # (20.0, -20.0)
from pyproj import CRS
CRS("EPSG:32632")
<Projected CRS: EPSG:32632>
Name: WGS 84 / UTM zone 32N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.
- bounds: (6.0, 0.0, 12.0, 84.0)
Coordinate Operation:
- name: UTM zone 32N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

The target CRS is EPSG:4326 which is a geographic CRS with units of degrees.

CRS("EPSG:4326")
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

As such, you are forcing the output resolution of the raster to 20 degrees. This is why you are getting a 1 x 1 raster.

The default projection seems fine:

>>> rds.rio.reproject(4326)
<xarray.DataArray (band: 1, y: 723, x: 951)>
array([[[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]]], dtype=uint8)
Coordinates:
  * x            (x) float64 12.44 12.44 12.44 12.44 ... 12.64 12.64 12.64 12.64
  * y            (y) float64 44.12 44.12 44.12 44.12 ... 43.96 43.96 43.96 43.96
  * band         (band) int64 1
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    _FillValue:    0