Is there anyway to output to grib2 after performing map algebra? For example, If you are performing a calculation between two grib2 bands, which rasterio does quite well, are you able to transform the numpy array into a grib2 file? I have been able to do this with geotiff but not grib2.

Ok, so it appears the issue was that I was attempting to do a windowed read outside the bounds of the geotiff (for TMS edge tiles) and `src.read()` was only returning the data inside the bounds since I wasn't using the `boundless=True` option. However, the next problem I am running into is that I'd like to use `fill_value=numpy.nan` to fill in the edge data with NaN values. However, my data is UInt8. Is there any way to specify an option to cast the data read by the call to `read()` into float32 so that I can support this?

Short version:

The following causes errors, in rasterio 1.0.15,  but is fine in rasterio 1.0.9. 

from rasterio.crs import CRS
print(CRS())

The following are the only tests in rasterio to check blank CRSs.

def test_null_crs_equality():
    assert CRS() == CRS()
    a = CRS()
    assert a == a
    assert not a != a


def test_null_and_valid_crs_equality():
    assert (CRS() == CRS(init='epsg:4326')) is False

I feel like print(CRS()) should not cause errors - CRS().__str__ should just return an empty string, for example.

Longer version:

We work with pre-orthorectified sat imagery, which have blank geotransforms (but may have GCPs set), eg,

gdalinfo /data/16APR27183709-P1BS-502376654080_01_P009.tif 
Driver: GTiff/GeoTIFF
Files: /data/16APR27183709-P1BS-502376654080_01_P009.tif
Size is 3201, 3281
Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
...

This causes problems on newer rasterio versions (presumably due to the switch to WKT representation for CRS). Example:

import rasterio
with rasterio.open('/data/16APR27183709-P1BS-502376654080_01_P009.tif') as src:
    print("CRS:", src.crs)

On rasterio 1.0.9, this is fine, src.crs evaluates to None:

python ~/test.py
/home/nick/venvs/ras9/lib/python3.6/site-packages/rasterio-1.0.9-py3.6-linux-x86_64.egg/rasterio/__init__.py:217: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: None

The NotGeoreferencedWarning is fine, we explicitly ignore it. The point is that src.crs is None.
On 1.0.15, however,

python ~/test.py
/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/__init__.py:216: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: Traceback (most recent call last):
  File "rasterio/_crs.pyx", line 274, in rasterio._crs._CRS.to_dict
  File "rasterio/_err.pyx", line 194, in rasterio._err.exc_wrap_ogrerr
rasterio._err.CPLE_BaseError: OGR Error code 7

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/nick/test.py", line 3, in <module>
    print("CRS:", src.crs)
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 245, in to_string
    return self.to_wkt() or self.to_proj4()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 106, in to_proj4
    return ' '.join(['+{}={}'.format(key, val) for key, val in self.data.items()])
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 170, in data
    self._data = self.to_dict()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 164, in to_dict
    return self._crs.to_dict()
  File "rasterio/_crs.pyx", line 277, in rasterio._crs._CRS.to_dict
rasterio.errors.CRSError: The WKT could not be parsed. OGR Error code 7

Interestingly, the GDAL Python API handles this in a different way:

from osgeo import gdal
ds = gdal.Open('/data/16APR27183709-P1BS-502376654080_01_P009.tif')
prj = ds.GetProjection()
print(prj)

which prints

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]

I guess GDAL assumes the GCP projection, in case of blank geotransform.

In any case, I think the src.crs is None behavior is good, since we can just look at src.gcps for the gcps.

I will add a workaround on our side, since

print(CRS().is_valid)

works, (that is, CRS().is_valid evaluates to False, without throwing exceptions.

Hi @Madhav, 
On linux it should work with `rasterio==1.0.15`