Problem with reproject epsg:3411 to epsg:4326

[trondkr]

Hi,

Thanks for a great Python package! I have a dataset I want to convert to a regular cartesian coordinate system. The files contain sea-ice concentration from the NSIDC whom provides the data as netCDF files. However, the projection of the data is in EPSG:3411. When converting these data using rioxarray reproject I run into a few issues and fail to create a valid file and I hope someone can help me get this correct.

I read the files (two example files can be found here:

                def preprocess(ds):
                    ds=ds.drop_vars(names=["melt_onset_day_cdr_seaice_conc_monthly",
                                       "nsidc_bt_seaice_conc_monthly","nsidc_nt_seaice_conc_monthly",
                                       "qa_of_cdr_seaice_conc_monthly","stdev_of_cdr_seaice_conc_monthly"])
                    ds = ds.swap_dims({"tdim": "time"})
                    units = "days since 1601-01-01 00:00:00"
                    dates = cftime.num2date(ds["time"].values, units,
                                        calendar="standard",
                                        only_use_cftime_datetimes=False)

                    if isinstance(dates[0], cftime.DatetimeGregorian):
                        d=dates[0]
                        dates=[datetime.datetime(d.year, d.month, d.day, d.hour, d.minute, d.second)]

                    return ds.assign({'time': (['time'], dates)})

                ds_in = xr.open_mfdataset("/CDR/seaice_conc_monthly_nh_*.nc",
                                            combine="by_coords",
                                            parallel=True,
                                            concat_dim="time",
                                            coords="minimal",
                                            data_vars="minimal",
                                            compat='override',
                                            preprocess=preprocess,
                                            decode_times=False,
                )
                ds_in = ds_in.rename({"cdr_seaice_conc_monthly": "siconc"})
                ds_in = xr.decode_cf(ds_in)
                ds_in = ds_in.rename({"xgrid": "x", "ygrid": "y"})
                ds_in = ds_in.reindex({"y":np.array(list(reversed(ds_in.y)))})
                ds_in = ds_in.transpose("time", "y", "x")

                ds_in = xr.decode_cf(ds_in)
                ds_in = ds_in.rio.write_crs("epsg:3411")

Then I rename xgridto x and ygridto have both dimensions and coordinates with the same name as expected by rioxarray. Next I convert the order of the y-dimension as data goes from north to south and I want south to north.

    ds_in = ds_in.rename({"xgrid": "x", "ygrid": "y"})
    ds_in = ds_in.reindex({"y":np.array(list(reversed(ds_in.y)))})

This creates a dataset that is supposedly an EPSG:3411 projection and makes sense:

Dimensions:     (time: 506, x: 304, y: 448)
Coordinates:
  * y           (y) float32 -5.338e+06 -5.312e+06 ... 5.812e+06 5.838e+06
    projection  int64 0
  * time        (time) datetime64[ns] 1978-11-01 1978-12-01 ... 2020-12-01
  * x           (x) float32 -3.838e+06 -3.812e+06 ... 3.712e+06 3.738e+06
Data variables:
    siconc      (time, y, x) float32 dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
...

Next I add the crs and reproject:

    ds_in_lonlat = ds_in.rio.reproject("epsg:4326")

This results in a new dataset where values, dimensions, and values are wrong:

ds_in_lonlat <xarray.Dataset>
Dimensions:     (time: 506, x: 1052, y: 75)
Coordinates:
  * x           (x) float64 -178.0 -177.7 -177.3 -177.0 ... 176.8 177.1 177.4
  * y           (y) float64 56.11 55.77 55.43 55.09 ... 32.09 31.76 31.42 31.08
  * time        (time) datetime64[ns] 1978-11-01 1978-12-01 ... 2020-12-01
    projection  int64 0
Data variables:
    siconc      (time, y, x) float32 2.54 2.54 2.54 2.54 ... nan nan nan nan
...

First, the x-dimension is now 1052 while y is 75 compared to earlier 304 and 448. I assume something in my approach is not correct and creates this erroneous result and I hope someone can help me find the error and get this conversion correct. I really appreciate it.

Cheers, Trond

[snowman2]

First, the x-dimension is now 1052 while y is 75 compared to earlier 304 and 448

When re-projecting, the shape is not usually the same. rasterio/GDAL automatically choose what it believes are the best dimensions in the new projection. If you want to preserve the same shape, you can pass it into the rio.reproject method.

[trondkr]

@snowman2 Thanks for your help with this. I added shape to the reproject function assuming shape=(y, x).

ds_in.rio.reproject("epsg:4326", shape=(448,304))

This provides a correctly sized dataset but the results look completely wrong. I have a feeling that one of the dimensions is rotated wrong, but can't see anything obviously wrong. If anyone has a ton of time and opportunity to take a look at the example file I would really appreciate it.

the current reprojection looks like this:

But should look more like this:

Thanks again.

[snowman2]

@trondkr, what happens when you do ds_in.rio.to_raster("file.tif") and then re-project with GDAL later?

[ppokhrel1]

@trondkr Did you solve it? I was running into the same problem.

Any help would be greatly appreciated.

Thanks for your time.

Pujan