issue_comments

I am not certain about my availability at the moment. If I get free time, I would be happy to do so. However, if someone else gets started on it first, that would be great as well.

Thanks @Illviljan :+1:

mypy should have caught this a while ago

Type hints are likely incomplete/inexact in some places in rioxarray code.

We have other plans to eventually make pandas an optional import too

That makes not using a data frame make a lot of sense.

I think show_engines provides a simple user friendly interface that is valuable for users in it's current state.

Looks like something already exists: python from xarray.backends.plugins import list_engines It would be nice to add the description & documentation properties to xarray.backends.common.BackendEntrypoint that can be implemented by external libraries (and internally) and a more official list_engines function with the above metadata could be added.

That's nice @dcherian :+1:

Example using nearest & tolerance with reindex & sel when dims don't match based on the example in sel:

```python import numpy import xarray

da = xarray.DataArray( numpy.arange(25).reshape(5, 5), coords={"x": numpy.arange(5), "y": numpy.arange(5)}, dims=("x", "y"), ) tgt_x = numpy.linspace(0, 4, num=5) + 0.5 tgt_y = numpy.linspace(0, 4, num=5) + 0.5 da = da.reindex( x=tgt_x, y=tgt_y, method="nearest", tolerance=0.2, fill_value=numpy.nan ).sel( x=xarray.DataArray(tgt_x, dims="points"), y=xarray.DataArray(tgt_y, dims="points"), ) Output: <xarray.DataArray (points: 5)> array([nan, nan, nan, nan, nan]) Coordinates: x (points) float64 0.5 1.5 2.5 3.5 4.5 y (points) float64 0.5 1.5 2.5 3.5 4.5 Dimensions without coordinates: points `` Side note: I don't think it makes sense to addfill_valuetoselas it would require adding new coordinates that didn't exist previously. Callingreindex` first makes that more clear in my opinion.

GDAL doesn't actually need the CRS stored with the _CRS convention as documented in their ZARR driver specification?

GDAL may need to update the ZARR driver. The NetCDF driver does support CF. It would be good for them to be consistent.

Top reasons for using CF over a CRS attribute:

• Supported by geospatial & netCDF software (GDAL/Esri/Unidata)
• Attributes are easily lost ref
• Encourage more usage of standard mechanisms of CRS storage

This document is also a useful reference for storing CRS in xarray: https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html

Instead, it recommends to use CF conventions for encoding CRS. This is more compatible with NetCDF, but won't be parsed correctly by GDAL.

GDAL does support the CF conventions for storing the CRS.

https://gdal.org/drivers/raster/netcdf.html#georeference

"The driver first tries to follow the CF-1 Convention from UNIDATA looking for the Metadata named “grid_mapping”."

https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html#xarray

Are you able to use decode_coords="all" when opening?

https://github.com/jupyter/nbconvert/releases/tag/6.4.3

No worries. I will try to address the conflicts when I get a free moment. Though, I wouldn't complain if someone else with a sooner free moment finished this one off before I do :-).

Sounds like that may need changes in rioxarray. Mind opening up an issue: https://github.com/corteva/rioxarray/issues ?

can we make it?

I am not seeing an option to do that: https://pre-commit.ci/

Is mypy being run in the CI now?

Yes, it runs on the CI

I didn't know we could pass --all like that.

I think it is equivalent to --all-files

This works for rioxarray. Thanks :+1:

I see the same issue reported here: https://github.com/pydata/xarray/issues/5944#issuecomment-962958896

See: https://github.com/corteva/rioxarray/runs/4140632105 (Note: updated CI build link)

When a fix is added for this, feel free to ping me to test the fix on rioxarray.

Thanks for the review 👍

Updated this so it is deprecated in 0.19.1 and added an entry to whats-new.rst.

Fix for doc build issues in #5805

xr.open_rasterio() is being used earlier in the pipeline

What issues do you have if you change to use rioxarray.open_rasterio?

This appears to not be an xarray issue since the code uses open_dataset and the rio.crs accessor. The rio.crs accessor is from rioxarray. From what I can tell, the issue has to do with the data is stored in the netCDF file: https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html

My guess is that there is an attribute named crs with the contents +init=EPSG:XXXX. If that is the case, I recommend using rioxarray to do:

python data.rio.write_crs("ESPG:XXXX", inplace=True)

Hopefully this is helpful.

For reference: - https://corteva.github.io/rioxarray/stable/examples/COG.html - https://github.com/corteva/rioxarray/pull/48

Also see: https://github.com/corteva/rioxarray/pull/359

Not sure what we want to do about these attributes in the long term.

I like the idea of using them the same way xarray does for metadata. The next step would be to get pandas to add .encoding and then they could do mask_and_scale on their read_* methods :).

Yes, the alignment stuff in there should do the trick.

I think it would be nice for xarray to have a "close enough" option where a certain number of decimal places is enough to consider them equal. Is there an issue somewhere where this is under consideration or is there an option already in existence that I missed?

Also related: https://github.com/corteva/rioxarray/discussions/332

@alexamici added the entrypoint here; https://github.com/corteva/rioxarray/pull/281. It seems xarray 0.18+ is required for it to work. Any timeline for the 0.18 release?

Related: https://github.com/corteva/rioxarray/issues/191

Similarly at https://github.com/corteva/rioxarray/issues/298#issuecomment-820557568

python np.abs(resampled.x.values - da2.x.values).max(), np.abs(resampled.y.values - da2.y.values).max() (4.000014541816199e-09, 6.666663665555461e-09)

For the original post for this issue ref: np.abs(raster_in.x.values - out_grid.x.values).max(), np.abs(raster_in.y.values - out_grid.y.values).max() (1.7763568394002505e-15, 7.105427357601002e-15)

Looks like there is a tiny difference in decimal precision.

Similar issue here: https://github.com/corteva/rioxarray/issues/298

Reproducible example in the issue ^^

Thanks for mentioning that. For rioxarray, the fix is: python transform = da.rio.transform()

For reference for how rioxarray does things: https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html

```

import xarray, rioxarray xda = xarray.DataArray(1) xda.rio.write_crs(4326, inplace=True) <xarray.DataArray ()> array(1) Coordinates: spatial_ref int64 0 Attributes: grid_mapping: spatial_ref xda.spatial_ref <xarray.DataArray 'spatial_ref' ()> array(0) Coordinates: spatial_ref int64 0 Attributes: crs_wkt: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["... semi_major_axis: 6378137.0 semi_minor_axis: 6356752.314245179 inverse_flattening: 298.257223563 reference_ellipsoid_name: WGS 84 longitude_of_prime_meridian: 0.0 prime_meridian_name: Greenwich geographic_crs_name: WGS 84 grid_mapping_name: latitude_longitude spatial_ref: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["... ```

Looks like entrypoint support was added in #4577. What is the timeline for the next release?

and rioxarray agrees

Is rioxarray agreeing to allow the addition of entrypoint support? If so, :+1:

You may be interested in rioxarray: https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html

Reference for writing from xarray to GeoTIFF or any other GDAL supported format: https://corteva.github.io/rioxarray/stable/examples/convert_to_raster.html

In rioxarray we updated open_rasterio so it opens a tif file like a DataArray and netcdf-like files as a dataset. It also loads in the metadata/tags as well:

NetCDF file as Dataset

python import rioxarray xds = rioxarray.open_rasterio("PLANET_SCOPE_3D.nc") <xarray.Dataset> Dimensions: (time: 2, x: 10, y: 10) Coordinates: * y (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06 * x (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05 * time (time) object 2016-12-19 10:27:29 2016-12-29 12:52:41.659696 spatial_ref int64 0 Data variables: blue (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368 green (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41 Attributes: coordinates: spatial_ref

GeoTiff as DataArray

python rds = rioxarray.open_rasterio("test_albedo.tif") <xarray.DataArray (band: 1, y: 936, x: 745)> [697320 values with dtype=int8] Coordinates: * band (band) int64 1 * y (y) float64 4.225e+06 4.225e+06 ... 4.216e+06 4.216e+06 * x (x) float64 6.168e+05 6.168e+05 ... 6.242e+05 6.242e+05 spatial_ref int64 0 Attributes: transform: (10.0, 0.0, 616800.0, 0.0, -10.0, 4225260.0) _FillValue: -1.0 scale_factor: 1.0 add_offset: 0.0 long_name: Albedo dry - the estimated ratio of the incident short-wav... grid_mapping: spatial_ref

So, it could be modified to use as an engine in its current state.

go for it!

Sounds good, will do. Thanks for the patch :+1:

I don't know much about rioxarray

The usage examples are a good place to go to learn more.

Is the intent to someday merge that into xarray or will the geospatial features mean it will remain separate?

It hasn't really been discussed/decided upon.

@cheginit, thanks for being willing. The part to modify is here: https://github.com/corteva/rioxarray/blob/6587726a019f46ec68ca94d6cea866abea956a2a/rioxarray/_io.py#L657

Would you like to add this update to rioxarray as well? If not, mind if I port it over?

Would you like to add this update to rioxarray as well? If not, mind if I port it over?

Would you like to add this update to rioxarray as well? If not, mind if I port it over?

Thanks for the fixes 👍

We might want to think about adding rioxarray as documentation build dependency

That brings along with it the fun of GDAL as a dependency.

@dcherian Maybe a link to these examples in that section then? https://corteva.github.io/rioxarray/stable/examples/examples.html

Adding an example on doing these manipulations would also be very welcome.

What types of manipulations are you thinking about? Reprojection? Export to geotiff? Clip? Reading in files?

Going to close this as it is over a year old.

For reference, WKT support is in rioxarray. https://corteva.github.io/rioxarray/stable/examples/crs_management.html

I'm wondering if an xarray accessor could work around some of the current limitations you're seeing.

I think this sounds like a good idea.

For reference, here is how CRS information is handled in rioxarray: CRS management docs.

@jsignell, the masked kwarg is in rioxarray if you need it: https://corteva.github.io/rioxarray/html/examples/clip_geom.html

Nice!

Put another way: why don't we put all the logic in rioxarray and make rioxarray an optional dependency of xarray to open rio files?

That is an option. All of the logic has already been moved over.

I'm curious, are there features in rioxarray that could be pushed upstream?

Depends on what the xarray maintainers would like to add. I would definitely like to see the open_rasterio function in rioxarray absorbed back into xarray. Things have just been really slow moving with the other PRs/issues.

I've abandoned this PR. If anyone has time to pick it up, that would be welcome.

I appreciate you staring this! Based on this PR, I added the feature into rioxarray here: https://github.com/corteva/rioxarray/pull/31 (released in version 0.0.9). (Example usage can be seen here.)

How to create tiled TIFF files in the tests

One thing I would like to note is that the automatic chunking would be useful if the raster is tiled or not. I tested out a raster that was not tiled, but it still had chunks. This is due to the raster being written in stripes. So, I would recommend removing the restriction to only tiled rasters.

Also, to create a tiled raster: ```python import rasterio import numpy from affine import Affine

with rasterio.open( "tiled.tif", "w", driver="GTiff", count=2, width=1024, height=1024, crs="+init=epsg:4326", transform=Affine(0.0083333333, 0.0, -180.00416666665, 0.0, -0.0083333333, 75.00416666665), dtype=rasterio.float32, tiled=True, blockxsize=512, blockysize=512, ) as rds: rds.write((numpy.random.rand(2, 1024, 1024)*10).astype(numpy.float32)) ```

Looks like they have this option in the tests:

python open_kwargs=dict( tiled=True, blockxsize=512, blockysize=512 ) with create_tmp_geotiff(nx=1024, ny=1024, nz=2, open_kwargs=open_kwargs) as (tmp_file, expected): ....

The right way to merge different dtypes and block shapes in the tiff file. Currently I'm assuming that they're uniform

I have yet to run into a raster that varies dtypes and block shapes across bands. Most of the time, they are single band rasters. And if they are not, they have had the same dtype and block shape. So, I think your assumption is a good one for most use cases.

Also, only a single dtype is allowed currently: https://github.com/pydata/xarray/blob/1d7bcbdc75b6d556c04e2c7d7a042e4379e15303/xarray/backends/rasterio_.py#L39-L40

A new project called rioxarray has a to_raster method with the default driver of GTiff.

You can use it like so:

``` import rioxarray import xarray

xds = xarray.open_rasterio("myfile.tif") wgs84_xds = xds.rio.reproject("EPSG:4326") wgs84_xds.rio.to_raster("myfile_wgs84.tif") `` It currently only supports 2d/3dDataArrays. So, you would have to iterate over your variables to export each one to a raster if you have aDataset`.

@fmaussion just a ping for what are your thoughts are for next steps on this?

Good question, I was thinking adding tests.

I think it would be good to write the dataset to a netcdf, and read it back in with xarray.open_dataset to ensure everything encodes/decides properly.

I definitely agree 💯 with you on avoiding redundancy. I actually had the same thought. However, part of the reason for using crs_wkt is motivated by CF 1.7 convention. I figured it was a small change that sacrifices redundancy for consistency across the community. The crs_proj4 name was just an attempt to be consistent with your original comment and the crs_wkt. But, I think renaming it to proj4 to remove redundancy is a good choice as well.

Tag: @fmaussion @djhoese

Since the main discussion is resolved, I will close this issue to prevent spam. I will open another issue for implementation details later.

Sounds like a good path forward to me.

Since the crs attribute will be deprecated, thoughts on replacing it with a crs coordinate and the PROJ.4 and WKT representations in the coordinate attributes?

Something like: <xarray.DataArray 'crs' ()> array(0) Coordinates: crs int64 0 Attributes: crs_wkt: PROJCS["UTM Zone 15, Northern Hemisphere",GEOGCS["WGS 84",D... crs_proj4: +proj=...

I think the solution would be to add a masked argument and property to the class RasterioArrayWrapper here.

And pass in the (masked=masked) argument here

That would pass riods.read( ..., masked=self.masked) when reading data from rasterio here

@jsignell, this looks like it would work - I like the implementation. Some things I would recommend to add is checks for when the nodata value varies across bands (maybe skip the mask and raise a warning) and when the nodata value is already NaN. Also, you should probably set the nodatavals to NaN if masked to resemble the behavior of xarray elsewhere.

I think @fmaussion would be a good reviewer for this PR.

It would definitely be a fantastic feature to have. Some discussion happened in #2308 to discuss standardizing methods for storing data that would definitely assist with this, however it did not seem like a priority at the time. Definitely feel free to share your opinion and use cases.

Feel free to re-open if this becomes something of interest. Thanks!

I was thinking just a "simple" library with a CRS object with utilities to create the object from WKT, PROJ.4, CF projection parameters, etc... And could take the object and convert it to any of the other formats. I wouldn't add any transform/resampling code and instead be a codebase that can be used by libraries that do those things. This way it would be useful to many different projects (geoxarray, pyresample, cartopy, metpy, ...). Thoughts?

Sounds good 👍. I see the CRS code potentially being useful outside of geoxarray. What are your thoughts on moving the CRS specific code into its own package?

Lots of good thoughts there.

I think a lot depends on who you plan on having for a user base. I like rasterio.crs.CRS, but it does require GDAL. I know GDAL is a heavy dependency (solves runtime problems and it a source of installation problems), but if your users are already using it, then it isn't too big of an ask. If not, you could look into something like pycrs. But, it looks like it hasn't been touched for a while (never a good sign). Or, there are other options, (requiring a more work & maintenance) such as re-creating the rasterio.crs.CRS object and using cython with a copy of only the spatial reference code from GDAL in the repo (might be better to create your own repo/package if you head in this direction). Or maybe someone could ask really nicely for the GDAL maintainers to consider making the spatial reference code a package on it's own (not likely, but would be really nice).

My preference to have the CRS object something created/retrieved by the accessor based on information in the file. If it is not, users will have to remove the CRS object when using to_netcdf() as the object will not be serializable and to_netcdf() will throw errors. I think allowing the user to set the CRS on the file (by adding the crs to .coords and grid_mapping to data_vars) by passing in a valid projection string (WKT, proj.4) geo.set_crs() would be a good idea. This way, when the user calls, to_netcdf() it will happily do what they want. In addition, the CRS will still be available later with the accessor when loading in with xr.open_dataset.

@fmaussion, I guess the main reason for the proposed change is to be able to read in the raster using open_rasterio() and be able to generate a netCDF using to_netcdf() that would be compatible with the rasterio model when using rasterio.open().

Also, since: "xarray.Dataset is an in-memory representation of a netCDF file", to me it makes more sense to store the data in a standard netCDF location that is consistent regardless of the backend. In doing so, other parts of xarray and libraries build off of xarray won't be required to find the same data in multiple locations.

Thanks for taking the time to consider this idea. I have used it quite a bit and it has made the datasets compatible across the various GIS tools I use (rasterio, GDAL, QGIS). It has definitely made my life easier and so I wanted to share what I have learned. But, it is just an idea. I really appreciate the work of the xarray developers and whatever y'all decide is fine with me.

The example I gave was just demonstrating that the dimension is not required for the crs coordinate variable.

I agree with the functionality that would support standardizing the crs as a coordinate variable on load for the case you specified. That way, the read in file will always behave the same and writing to netCDF it would be output correctly.

This all sounds like it is heading in a good direction. 👍

It is not in the dimension, it is the coordinate attribute in the variable. That is handled automatically by xarray when writing to_netcdf if you add the crs to the coords.

From the ncdump: ``` dimensions: x = 65 ; y = 31 ; variables: double x(x) ; x:_FillValue = NaN ; x:long_name = "x coordinate of projection" ; x:standard_name = "projection_x_coordinate" ; x:units = "m" ; double y(y) ; y:_FillValue = NaN ; y:long_name = "y coordinate of projection" ; y:standard_name = "projection_y_coordinate" ; y:units = "m" ; int64 time ; time:units = "seconds since 2015-04-03T17:55:19" ; time:calendar = "proleptic_gregorian" ; int64 spatial_ref ; spatial_ref:spatial_ref = "PROJCS[\"UTM Zone 15, Northern Hemisphere\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-93],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]" ; double ndvi(y, x) ; ndvi:_FillValue = NaN ; ndvi:grid_mapping = "spatial_ref" ; ndvi:coordinates = "spatial_ref time" ;

// global attributes: :creation_date = "2018-04-11 13:14:55.401183" ; ```

It would definitely be a good idea to ensure that the crs variable is a coordinate, so I agree that having support for that would be a good idea.

That is interesting, I am definitely not an expert with non-uniform datasets. From the satellite datasets I have used, the 2D latitude and longitude coordinates are stored in the datasets and are not super useful. I usually have to use other ways to recreate the grid coordinates in the original projection (ex. SMAP uses the EASE Grid 2.0 but it stores the latitude/longitude of the points in the file) or reproject & flatten the coordinates. I have had to do this with weather data and made an xarray extension pangaea to handle it. So, that is what I was referring to when I misunderstood your question.

For your example of adding a crs attribute, ...

The files I have created have the crs coordinate variable inside the netCDF file already and it is always there when I load it in with xarray.open_dataset(). The method set_crs() could be used to add the crs coordinate variable and grid_mapping attributes to the dataset in the proper way so that it would be there on xarray.open_dataset() after dumping it to the file with to_netcdf().

The CF stuff is supported by rasterio, GDAL, QGIS and that is why I like it. If there is another way that is as well supported, I am not opposed to it.

In your example of methods is to_projection a remapping/resampling operation? If not, how does it differ from set_crs?

The to_projection() method would be a reproject/resampling operation.

That is strange. Was the variable put in data_vars and not in coords?

In your own projects and use of raster-like data, do you ever deal with non-uniform/non-projected data?

I have dealt with non-uniform data in the geographic projection. I have found it easiest to deal with it if you can determine the original projection and project the coordinates back to that projection so it is uniform. But, I am by no means an expert in this arena. Most if the time I work "normal" data.

How do you prefer to handle/store individual lon/lat values for each pixel?

rasterio/GDAL/QGIS all seem to use the centroid.

Also it looks like xarray would have to be updated to add the "crs" coordinate since currently it is not considered a coordinate variable. So a new library may need to have custom to_netcdf/open_dataset methods, right?

Actually, it is not difficult to add as it stands: python ds.coords['crs'] = 0 ds.coords['crs'].attrs = dict(spatial_ref="PROJCS["UTM Zone 15, Northern Hemisphere",GEOGCS["WGS 84",D...") But, if a crs does not already exist on the dataset, I guess that a function that adds the crs properly would be useful so it can also add the grid_mapping to all of the variables.

Example: python ds.geo.set_crs("+init=epsg:4326")

I think that minor modifications will be needed once the crs is set properly on the xarray dataset. Because after that, the to_netcdf() will automatically produce georeferenced datasets.

I could see the first pass of the extension/library simply performing: 1. ds.geo.crs 2. ds.geo.set_crs() 3. ds.geo.to_projection()

I don't think so. The crs should already exist on the dataset beforehand and I think it may be a bit too much to expect it to add it to the dataset in that function as there are some that don't have it and not all xarray datasets are geospatially located.

Here is an example of how it would look on a dataset: <xarray.Dataset> Dimensions: (x: 65, y: 31) Coordinates: * x (x) float64 ... * y (y) float64 ... time datetime64[ns] ... crs int64 ... Data variables: ndvi (y, x) float64 ... Attributes: Here is how the crs or spatial_ref coodinate variable would look: <xarray.DataArray 'crs' ()> array(0) Coordinates: time datetime64[ns] ... crs int64 0 Attributes: spatial_ref: PROJCS["UTM Zone 15, Northern Hemisphere",GEOGCS["WGS 84",D... And here is how it would look on the variables: <xarray.DataArray 'ndvi' (y: 31, x: 65)> array([[ ...]]) Coordinates: * x (x) float64 ... * y (y) float64 ... time datetime64[ns] ... crs int64 0 Attributes: grid_mapping: crs

@djhoese Whether or not we use the CF convention is not what I am concerned about. What I think would benefit the most people is with the file format to be able to do to_netcdf() with the file and be able to have it read in with standard GIS tools such as rasterio, GDAL, and QGIS. With this schema, this is possible.

Another benefit is that it keeps the crs or spatial_ref with it when you do your operations as it is a coordinate of the variable.

Also, as a side note if you use the center pixel coordinates, then GDAL, rasterio, and QGIS are able to read in the file and determine it's affine/transform without a problem.

For the new library, if you have a crs method attached to it, it isn't too difficult to convert it to whatever format you need it to be in. With the rasterio.crs.CRS.from_string you can "Make a CRS from an EPSG, PROJ, or WKT string" and you can also get back the EPSG, PROJ, or WKT string with a simple method call.

For example, using the recommended method to extend xarray, you could add a crs property:

python from rasterio.crs import CRS ........ @property def crs(self): """:obj:`rasterio.crs.CRS`: Projection from `xarray.DataArray` """ if self._crs is not None: return self._crs try: # look in grid_mapping self._crs = CRS.from_string(self._obj.coords[self._obj.grid_mapping].spatial_ref) except AttributeError: raise ValueError("Spatial reference not found.") return self._crs

And if you call your extension geo, all you would need to get the CRS would be: python ds.geo.crs To get proj.4 string: python ds.geo.crs.to_string() To get WKT string: python ds.geo.crs.wkt To get EPSG code: python ds.geo.crs.to_epsg()

I am really excited about this discussion. I know of other libraries that have done the same thing and have written internal libraries myself.

If possible, I would hope that we could follow the CF convention on this as it makes the output netCDF file compatible with QGIS, GDAL, and rasterio when written using to_netcdf()

To do so, you add the crs coordinate to the dataset/dataarray. int crs; And then, you add the spatial_ref attribute to the crs which is a crs WKT string.

Next, you add the grid_mapping attribute to all associated variables that contains crs as the grid_mapping.

See an example here.

After that, you could store all kinds of information inside the crs variable such as the proj4 string, the affine, etc.

Indeed, it would be nice to refactor the tests to use pytest functions so we can use fixtures and the pytest.mark.parametrize to reduce duplication.

I second this request. I have had this same need before where I don't always know what variables I don't want, but I do know the variables that I do want.

Would it be possible to put the coordinate generation into a separate function? It would be nice to pass in the affine, width, height, and is_2d as args and get back the spatial coordinates.

One thought I had was that there could be a similar parameter to mask_and_scale added that could be just called mask so open_rasterio behaves similarly to the open_dataset api. The mask parameter can allow the user to decide if they want to convert to float to get a mask or leave it alone if they want to preserve the integer data.

I would vote for adding the _FillValue attribute as it conforms with the CF conventions and it will make the dataset work better with other software looking for this value if you write it out with to_netcdf().

Also, on a similar vein, it would be nice to break out the bands into separate variables if it is a multi-band raster. This way, each band can store it's own nodata value and metadata separately in the dataset.

@fmaussion, thanks again for your time!

Which is what I think you do in salem.

Oh, and to make the time slice work for WRF, I had to make the time variables consistent: python xds.rename( { time_dim: 'time', time_var: 'time', }, inplace=True )

OK. That is good to know. In this case, I don't know the name if the dimension beforehand, so I thought it might be preferable to do:

python data = data.loc[{self.lsm_time_dim: [pd.to_datetime(time_step)]}] instead of: python data = data.sel(**{self.lsm_time_dim: [pd.to_datetime(time_step)]}) What are your thoughts?

OK. I feel silly now. I thought I should have been using loc[] Thanks @fmaussion for that.

I am having the same issue with a completely different grid: <xarray.DataArray 'sp' (time: 72, latitude: 5, longitude: 4)> dask.array<getitem, shape=(72, 5, 4), dtype=float64, chunksize=(24, 5, 4)> Coordinates: * longitude (longitude) float32 248.0 248.25 248.5 248.75 * latitude (latitude) float32 40.75 40.5 40.25 40.0 39.75 * time (time) datetime64[ns] 2016-01-02 2016-01-02T01:00:00 ... Attributes: units: Pa long_name: Surface pressure standard_name: surface_air_pressure Error: ``` self = DatetimeIndex(['2016-01-02 00:00:00', '2016-01-02 01:00:00', '2016-01-02 02:00:00', '2016-01-02 03:00:0...:00:00', '2016-01-04 22:00:00', '2016-01-04 23:00:00'], dtype='datetime64[ns]', freq=None) key = array(Timestamp('2016-01-02 00:00:00'), dtype=object)

def __getitem__(self, key):
    """
        This getitem defers to the underlying array, which by-definition can
        only handle list-likes, slices, and integer scalars
        """

    is_int = is_integer(key)
    if is_scalar(key) and not is_int:
        raise ValueError

    getitem = self._data.__getitem__
    if is_int:
        val = getitem(key)
        return self._box_func(val)
    else:
        if com.is_bool_indexer(key):
            key = np.asarray(key)
            if key.all():
                key = slice(0, None, None)
            else:
                key = lib.maybe_booleans_to_slice(key.view(np.uint8))

        attribs = self._get_attributes_dict()

        is_period = isinstance(self, ABCPeriodIndex)
        if is_period:
            freq = self.freq
        else:
            freq = None
            if isinstance(key, slice):
                if self.freq is not None and key.step is not None:
                    freq = key.step * self.freq
                else:
                    freq = self.freq

        attribs['freq'] = freq

      result = getitem(key)

E IndexError: arrays used as indices must be of integer (or boolean) type

../../../tethys/miniconda/envs/gssha/lib/python3.6/site-packages/pandas/core/indexes/datetimelike.py:296: IndexError ```

I just did. Thanks!

Strange. But I can't seem to re-produce the issue. Maybe it was on a Windows machine or maybe it is fixed now.

@shoyer, thanks for looking into it. I am resampling from 3hr data to 1hr data.

python resampled_ds = ds.resample('1H', dim='time', keep_attrs=True)

I am using it here: https://github.com/CI-WATER/gsshapy/blob/f4e5cb13c1d528021e1953859b712553a4162311/gsshapy/grid/grid_to_gssha.py#L789-L844

I ran into the issue there and had to add code to make sure the coordinates were copied.

Thanks!

I did a $ ncdump -v time file.nc And I got: ``` ... float time(time) ; time:units = "minutes since 2011-03-05 03:00:00" ; time:long_name = "time" ; time:time_increment = "10800" ; time:begin_date = "20110305" ; time:begin_time = "030000" ; ... data:

time = _ ; } ```

So, it looks like a bad file ...