0

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?

arctic_climate_science
  • 57
  • 5
  • 1
    best to reproject to a target rather than heuristics for a given crs (this case should be ok, but the reverse has no bounds in stereographic when software tries to match the resolution rather than the dimension), set the extent and dimension explicitly to the lat and lon range you want - set the lat range to 0,90, say – mdsumner Jan 05 '23 at 02:46
  • 1
    fwiw I tried this with gdalwarp directly and it makes good choices, you get lat [0,90] - so maybe rioxarray can fall back to what the warper does here, seems like worth investigating for the project – mdsumner Jan 05 '23 at 02:58
  • using: `gdalwarp -te -180 0 180 90 -t_srs EPSG:4326 cmc_sdepth_dly_1998_v01.2.tif cmc_sdepth_dly_1998_v01.2_latlon.tif` looks like it may have worked! Thank you! – arctic_climate_science Jan 05 '23 at 17:42

0 Answers0