1

I am getting different resutls when plotting a np.array with matplotlib than when doing it in an external software - QGIS. Here is what I have:

    print(wiw[0].shape)
(1, 8433, 9254)

    print(wiw[0].min())
nan

    print(wiw[0].max())
nan

    print(wiw[0].dtype)
float64
    print(np.nanmin(wiw[0]))
1.0
    print(np.nanmax(wiw[0]))
1.0

This is how it looks like when I plot it with matplotlib:

fig, ax = plt.subplots(figsize = (20,20))
ax.imshow(wiw[0][0,:,:], cmap='Greys_r');

enter image description here

Next, I save it with rasterio using the following code:

with rasterio.open(output_path + 'WIW_TEST.tif', 'w', **ras_meta) as dst:
    dst.write(wiw[0])

Where ras_meta is:

{'driver': 'GTiff', 'dtype': 'float64', 'nodata': 0, 'width': 9254, 'height': 8433, 'count': 1, 'crs': CRS.from_epsg(32633), 'transform': Affine(20.0, 0.0, 300000.0,
       0.0, -20.0, 1600020.0), 'blockxsize': 640, 'blockysize': 640, 'tiled': True}

If I open this geotiff in QGIS, its appearence is differnet:

enter image description here

But surprisingly, if I open the saved geotiff again with rasterio and plot it with matplotlib, the visualization is correct. What is going on? Is it a dtype issue, a nodata problem?

4
  • 2
    QGIS and matplotlib both have to stretch/scale the data for display, they're just doing it slightly differently. You can tweak the stretch in QGIS layer symbology if you want.
    – user2856
    Sep 20 at 12:29
  • Could be but in this case my array is formed by np.nan and 1. There are not intermediate values. For reference, this array was created with a np.where and a condition. Where it was met, 1. If not, a np.nan was put
    – GCGM
    Sep 20 at 13:52
  • Why are you storing binary data as a float?
    – mikewatt
    Sep 21 at 0:35
  • 1
    If I am not wrong, in order to have nan, float is needed as dtype
    – GCGM
    Sep 21 at 6:16
1

I suspect this has to do with the styling when drawing on the QGIS canvas, not with any data differences. Have you tried changing your resampling options in the symbology panel of QGIS? For example, change from nearest neighbor to bilinear and tinker with the oversampling to see if that matches what you expect more closely.

You can also try plotting the histogram of the raster when you have it in QGIS to verify the underlying data is being read correctly. If the problem is purely how the raster looks in the QGIS canvas, and you can verify the underlying data is the same, I would not sweat it too much.

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.