I am attempting to reproject a CMC Snow Depth Map (in GeoTIFF format) from North Polar Stereo to Lat/Lon. The PROJ string of the file is: '+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs=True'. The data is stored here: https://drive.google.com/file/d/1EhzWiOMF2YEYpDOPpmfH58tJMQOxNdIM/view?usp=share_link
Here is a reproduction of my code block:
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.crs import CRS
from rioxarray.crs import crs_from_user_input
import rioxarray as rxr
#### Open CMC Data ####
yrs = np.arange(1980,2020,1)
CMC_dir = '/mnt/data/obs/CMC_snow_depth/'
CMC_fil = ''.join([CMC_dir+'cmc_sdepth_dly_1998_v01.2.tif'])
dates = pd.date_range(start='08-01-1998', end='12-31-1998', freq='D')
### Create XArray Dataset ###
CMC_dataset = rxr.open_rasterio(CMC_fil,decode_coords="all")
print(CMC_dataset)
#writing the crs, using proj4
crs_CMC = CRS.from_string('+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs=True')
CMC_dataset = CMC_dataset.rio.write_crs(crs_CMC)
### Reproject to Lat/Lon from Polar Stereographic ###
CMC_dataset_latlon = CMC_dataset.rio.reproject("EPSG:4326") #EPSG:4326 is the lat/lon projection
print(CMC_dataset_latlon)
The input data has the following structure:
<xarray.DataArray (band: 153, y: 706, x: 706)>
[76260708 values with dtype=float64]
Coordinates:
* band (band) int64 1 2 3 4 5 6 7 8 ... 147 148 149 150 151 152 153
* x (x) float64 -8.394e+06 -8.37e+06 ... 8.37e+06 8.394e+06
* y (y) float64 8.394e+06 8.37e+06 ... -8.37e+06 -8.394e+06
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
_FillValue: -1.7e+308
scale_factor: 1.0
add_offset: 0.0
However, when I go to reproject the file to Lat/Lon, the dataset becomes corrupted, and the lat/lon become strange:
<xarray.DataArray (band: 153, y: 54, x: 997)>
array([[[ 0.0e+000, 0.0e+000, 0.0e+000, ..., 0.0e+000, 0.0e+000,
0.0e+000],
[ 0.0e+000, 0.0e+000, 0.0e+000, ..., 0.0e+000, 0.0e+000,
0.0e+000],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., 0.0e+000, 0.0e+000,
-1.7e+308],
...,
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308]],
[[ 0.0e+000, 0.0e+000, 0.0e+000, ..., 0.0e+000, 0.0e+000,
0.0e+000],
[ 0.0e+000, 0.0e+000, 0.0e+000, ..., 0.0e+000, 0.0e+000,
0.0e+000],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., 0.0e+000, 0.0e+000,
-1.7e+308],
...
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308]],
[[ 0.0e+000, 0.0e+000, 0.0e+000, ..., 0.0e+000, 0.0e+000,
0.0e+000],
[ 0.0e+000, 0.0e+000, 0.0e+000, ..., 0.0e+000, 0.0e+000,
0.0e+000],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., 0.0e+000, 0.0e+000,
-1.7e+308],
...,
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308],
[-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
-1.7e+308]]])
Coordinates:
* x (x) float64 -179.8 -179.5 -179.1 -178.7 ... 179.1 179.4 179.8
* y (y) float64 19.3 18.94 18.57 18.21 ... 0.8805 0.5194 0.1583
* band (band) int64 1 2 3 4 5 6 7 8 ... 147 148 149 150 151 152 153
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0
_FillValue: -1.7e+308
Any idea what is happening here?