to_raster doesn't maintain CRS for xarray 0.19 #417
Comments
|
@JackLidge can you upgrade the version of rioxarray to the latest version and confirm whether this issue still persists? |
|
Yep, I upgraded to 0.7.1 and the same issue still occurred. Weird that conda was picking such an old version of rioxarray to install. |
|
What happens when you do this? wkt = 'PROJCS["WGS_1984_UTM_Zone_47N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",99.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'
ds.rio.write_crs(wkt, inplace=True)
ds = ds.rio.reproject('+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs', resolution=res)
ds = xr.where((ds == -9999) | (ds == 0), 15000, ds)
os.makedirs(os.path.dirname(tmp_fname), exist_ok=True)
ds.rio.to_raster(f'{tmp_fname}', dtype=np.int16, nodata=15000) |
|
Also, what does |
|
ds_out is just a DataArray containing the same data in ds. Print out of ds_out below: I think previously it flagged an error using to_raster on a dataset, so I converted a dataset to a dataarray when using to_raster. Probably an issue with previously using such an old version of rioxarray. No issues with writing out a raster just using ds, but still no CRS information being provided from gdalinfo. |
|
Are you able to provide the source raster? |
|
EDIT: File upload now below. This is the code I've been using to open the file to get it into an xarray reading for reading out: |
|
Side question: Is there a reason you aren't using ds = rioxarray.open_rasterio(in_fname)Your logic seems like the reason is related to #296. However, the final dataset doesn't match that, so maybe not. |
|
I wrote the script which does this processing a while ago - at the time this was the only method I could find to get everything to work in the way it does currently. I can't remember exactly why I used this specific combination of gdal, xarray and rioxarray to do everything. At the moment it works fine by just holding xarray at v0.16 but if this isn't can't be solved for whatever reason then I'll probably go back and rewrite it to see if it will work using |
|
It seems I may need that With osgeo.gdal: With rasterio: My environment: |
|
I am wondering if it is a file format where there are multiple files and some of them are missing in the zip you uploaded? |
|
There is an .hdr file which might be needed to download it, although I've never explicited checked whether it's needed for loading it or not. So I've re-uploaded the .bsq file with its corresponding .hdr file. I've also added the environment.yml file used in my conda environment. |
|
This works for me: import rasterio
import rioxarray
with rasterio.open("T47TQG_20210919T035541_TOA_rad_10m_atm.bsq") as rds, rasterio.vrt.WarpedVRT(rds, crs="EPSG:32648") as wds:
xds = rioxarray.open_rasterio(wds)
original_crs = xds.rio.crs
xds = xarray.where((xds == -9999) | (xds == 0), 15000, xds)
xds.rio.write_crs(original_crs, inplace=True)
xds.rio.write_nodata(15000, inplace=True)
xds.rio.to_raster('out.tif', dtype="int16")out_xds = rioxarray.open_rasterio("out.tif")
print(out_xds.rio.crs) |
|
Updated for: https://corteva.github.io/rioxarray/stable/getting_started/manage_information_loss.html import rasterio
import rioxarray
with rasterio.open("T47TQG_20210919T035541_TOA_rad_10m_atm.bsq") as rds, rasterio.vrt.WarpedVRT(rds, crs="EPSG:32648") as wds:
xds = rioxarray.open_rasterio(wds)
original_attrs = xds.attrs.copy()
original_encoding = xds.encoding.copy()
original_crs = xds.rio.crs
xds = xarray.where((xds == -9999) | (xds == 0), 15000, xds)
# preserve lost informaiton
xds.rio.update_attrs(original_attrs, inplace=True)
xds.rio.update_encoding(original_encoding, inplace=True)
xds.rio.write_crs(original_crs, inplace=True)
xds.rio.write_nodata(15000, inplace=True)
# write to file
xds.rio.to_raster('out_attrs.tif', dtype="int16") |
|
Note: The |
|
Yep, can confirm that works. Thanks very much for all your help, I'll look at reconfiguring this script to use that method. I haven't seen the vrt.WarpedVRT function before, is the CRS supplied there the input or output CRS? I wasn't fully sure based off scanning the rasterio docs for the function. |
|
Glad to hear that it worked. The crs kwarg is the destination CRS. Here are the API docs: https://rasterio.readthedocs.io/en/latest/api/rasterio.vrt.html?highlight=warpedvrt#rasterio.vrt.WarpedVRT |
|
Going to close as this issue seems to be resolved. |
snowman2
added
question
I've got some code which reads in data from .BSQ files using gdal, creates an xarray dataset and then reprojects and writes it out using rioxarray. This used to work fine but recently the CRS information hasn't been present in the written out datasets. When I downgrade the xarray version, this then works fine. I'm wondering if the way to assign a CRS to a dataset has changed (either in xarray or rioxarray I'm not sure) between version 0.16 and 0.19.
An example of the code being used is below (I can provide a full reproduceable code sample if required, but thought I'd check I'm not using some functions incorrectly first).
Code Sample, a copy-pastable example if possible
Problem description
This produces GeoTiffs with no CRS information in them. Running
gdalinfoon a dataset written out using the above method (xarray v0.19.0) results in the following output:Expected Output
Downgrading to xarray v.0.16.2:
Environment Information
For
xarray=0.16.2:For
xarray=0.19:Installation method
Conda install from environment.yml file, can provide full
conda listoutput orenvironment.ymlparameters if required, but both are quite long.The text was updated successfully, but these errors were encountered: