Unable to write raster from array using rasterio

arraygdalpythonrasterrasterio

I have two raster files, raster1 and raster2. I want to take raster1 and divide it by raster2 and produce an output raster for this calculation. I am using python and the rasterio package.

For info on raster1, I run: raster1.profile and receive:

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': None, 'width': 31232, 'height': 25630, 'count': 1, 'crs': CRS.from_wkt('LOCAL_CS["NAD83 / Conus Albers",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5070"]]'), 'transform': Affine(25.0, 0.0, 5983880.159924698,
       0.0, -25.0, 2476597.307965625), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'interleave': 'band'}

I then run raster2.profile and receive:

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': None, 'width': 31232, 'height': 25630, 'count': 1, 'crs': CRS.from_wkt('LOCAL_CS["NAD83 / Conus Albers",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5070"]]'), 'transform': Affine(25.0, 0.0, 5983880.159924698,
       0.0, -25.0, 2476597.307965625), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'interleave': 'band'}

Looks to be the exact same information.

I run the following to carry out my calculation:

raster1 = r"raster1.tif"
raster1 = gdal.Open(raster1)
rasterArray1 = raster1.ReadAsArray()

raster2 = r"raster2.tif"
raster2 = gdal.Open(raster2)
rasterArray2 = raster2.ReadAsArray()

new_raster = rasterArray1/rasterArray2
new_raster

And I see this output:

array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

(Not sure if all those nan values are an issue, but I will proceed)

I then want to convert this array into a raster .tif file that I can view in QGIS. I know that rasterio requires a raster to serve as the source for all the raster creation details, and so I use raster 2 for that and run:

src = rasterio.open('raster2.tif')
profile = src.profile
with rio.open('test.tif', 'w', **profile) as dst:
    dst.write(new_raster.astype(rio.uint8), 1)

This runs successfully and I see test.tif produced. Though when I go to look at this output raster in QGIS, I just see a black square with the legend showing the values going from 0 to 0. And then I see a QGIS error reading: No transform is available between Unknown CRS: ENGCRS["NAD83 / Conus Albers",EDATUM[""],CS[Cartes… and Unknown CRS: GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["WGS 8…. No coordinate operations are available between these two reference systems.

I am not sure what is going wrong, since it seems like my raster files are OK, and that I entered all the correct arguments necessary to divide arrays and then write the output array to a new raster. How can I fix my calculations and raster processing code so that I can properly create this new output raster?

Best Answer

You're writing a floating point array as an unsigned 8 bit integer. Try

src = rasterio.open('raster2.tif')
profile = src.profile
with rio.open('test.tif', 'w', **profile) as dst:
    dst.write(new_raster, 1) # No need to cast to int, as profile is already float32

or if that doesn't work, enforce the casting using .astype(rio.float32) instead of your original uint8

Related Question