2

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?

1 Answer 1

5

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

2
  • Thanks. I tried your first suggestion and it works, and I see my newly calculated raster displayed in QGIS. Though when I open this raster in QGIS I see an error message 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 if this is an issue or how I would address this. Jan 17 at 20:13
  • Check what the projection of the output raster is with e.g. gdalinfo & check it's the same as the source files.
    – Jose
    Jan 18 at 11:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.