I have a geotiff and a data array with the same shapes. The data array is auxiliary data of the geotiff. The x and y coordinates are the size of the geotiff coordinates, but the data array has no projection. I am trying to project the dataarray with the information from the geotiff.
My idea was:
- Make a new xarray data set from the above
- Reproject the new data set (to be able to plot it later on with a third data set on top)
The issues I have:
Using xarr.rio.reproject() is very smart, but my new data set needs a crs, which it do not have. I have tried:
#!/usr/bin/env python3
import numpy as np
import xarray as xr
import pyproj
import cartopy.crs as ccrs
import rioxarray as rio
from rasterio.warp import transform
xarr = xr.open_rasterio('my_geotif.tif')
ds = xr.open_dataset('my_dataarray.h5')
# Projection information
proj_tiff = xarr.crs
wgs84_epsg = 'epsg:4326'
proj_stereo = 'epsg:3995'
print(proj_tiff)
Original projection:
+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs=True
### Making the new data set with coordinates from tiff_file and data from dataarray
ds_new = xr.DataArray(ds.classification, coords={"y": xarr.y, "x": xarr.x}, dims=['y', 'x'])
ds_new.attrs = ds.attrs
# Adding CRS information
ds_new['crs'] = 0
ds_new['crs'].attrs = xarr.rio.crs
## closing used datasets
xarr.close()
ds.close()
## Reprojecting
ds_reproj = ds.rio.reproject(proj_stereo)
I get the error message:
File "/g5/procdata/skr/anaconda3/envs/pygmt/lib/python3.8/site-packages/rioxarray/rioxarray.py", line 1169, in reproject
raise MissingCRS(
rioxarray.exceptions.MissingCRS: CRS not found. Please set the CRS with 'set_crs()' or 'write_crs()'. Data variable: classification
I have tried 'set_crs() and write_crs(), but then I get: AttributeError: 'DataArray' object has no attribute
Any idea of what to do? It is probably simple, but I have tried so much by now, and nothing seams to work.
