question: save geotiff without modifying pixel value

[ybcheng]

Hi,

I'm using satpy to read, resample and save GOES-16 L1b radiance data. I've noticed the pixel values in the output geotiff file are re-mapped to the new bit-depth (e.g. 8bit or 16bit) instead of the original value. For instance, if the max pixel value is 200 in a scene, it will be scaled to 255 in a 8bit image or 65535 in a 16bit image. Is there any way to disable the auto scaling of pixel values? the values in the radiance data has physical meanings and I don't them to be changed.

Thank you,

[djhoese]

This isn't terribly well documented, but see the docstring for the geotiff writer here. In summary:

scn.save_datasets(writer=’geotiff’, dtype=np.float32, enhance=False)

[mraspaud]

@ybcheng Did you get it to work ?

[dtlloyd]

Hi,

I am piggy-backing on to this issue, because I have a similar problem. I would like to preserve the values in an image from dataset. At the moment when I use save_datasets or save_dataset, the image is changed so that the pixel values lie between 0 and 255.

Running the line of code djhoese recommended above I get a message:
"ValueError: cannot convert float NaN to integer"
I am using SatPy with a Sentinel 3 dataset. I get the same error on a dataset regardless of whether it has been resampled or not. Any ideas what might be causing this?

[djhoese]

@dtlloyd Usually in cases like this it is best to provide the full traceback for your error so we can see where the issue is coming from.

[dtlloyd]

@djhoese Thanks for the quick reply.
After running: newscn.save_datasets(writer='geotiff', dtype=np.float64, enhance=False)

The full traceback looks like:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-377-35894c5bb653> in <module>
     24 
     25 newscn.load([filename_S8])
---> 26 newscn.save_datasets(writer='geotiff', dtype=np.float64, enhance=False)
     27 #newscn.save_dataset(filename_S8)
     28 #newscn.show(filename_S8) # saves to

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/satpy/scene.py in save_datasets(self, writer, datasets, compute, **kwargs)
   1179                                "dimensions (eg. through resampling).")
   1180         writer, save_kwargs = load_writer(writer, ppp_config_dir=self.ppp_config_dir, **kwargs)
-> 1181         return writer.save_datasets(datasets, compute=compute, **save_kwargs)
   1182 
   1183     @classmethod

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/satpy/writers/__init__.py in save_datasets(self, datasets, compute, **kwargs)
    637         results = []
    638         for ds in datasets:
--> 639             results.append(self.save_dataset(ds, compute=False, **kwargs))
    640 
    641         if compute:

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/satpy/writers/__init__.py in save_dataset(self, dataset, filename, fill_value, overlay, decorate, compute, **kwargs)
    768         img = get_enhanced_image(dataset.squeeze(), enhance=self.enhancer, overlay=overlay,
    769                                  decorate=decorate, fill_value=fill_value)
--> 770         return self.save_image(img, filename=filename, compute=compute, fill_value=fill_value, **kwargs)
    771 
    772     def save_image(self, img, filename=None, compute=True, **kwargs):

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/satpy/writers/geotiff.py in save_image(self, img, filename, dtype, fill_value, floating_point, compute, **kwargs)
    246             # we can use the faster rasterio-based save
    247             return img.save(filename, fformat='tif', fill_value=fill_value,
--> 248                             dtype=dtype, compute=compute, **gdal_options)
    249         except ImportError:
    250             LOG.warning("Using legacy/slower geotiff save method, install "

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/trollimage/xrimage.py in save(self, filename, fformat, fill_value, compute, keep_palette, cmap, **format_kwargs)
    272         else:
    273             return self.pil_save(filename, fformat, fill_value,
--> 274                                  compute=compute, **format_kwargs)
    275 
    276     def rio_save(self, filename, fformat=None, fill_value=None,

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/trollimage/xrimage.py in pil_save(self, filename, fformat, fill_value, compute, **format_kwargs)
    391             format_kwargs['pnginfo'] = self._pngmeta()
    392 
--> 393         img = self.pil_image(fill_value, compute=False)
    394         delay = img.save(filename, fformat, **format_kwargs)
    395         if compute:

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/trollimage/xrimage.py in pil_image(self, fill_value, compute)
    658 
    659         """
--> 660         channels, mode = self.finalize(fill_value)
    661         res = channels.transpose('y', 'x', 'bands')
    662         img = dask.delayed(PILImage.fromarray)(np.squeeze(res.data), mode)

~/anaconda3/envs/No_Pydrive3pt6/lib/python3.6/site-packages/trollimage/xrimage.py in finalize(self, fill_value, dtype, keep_palette, cmap)
    640                     final_data = final_data.where(final_data != ifill, dtype(fill_value))
    641                 elif fill_value is not None:
--> 642                     final_data = final_data.fillna(dtype(fill_value))
    643 
    644         final_data = final_data.astype(dtype)

ValueError: cannot convert float NaN to integer

[djhoese]

This looks like it is trying to save to a PNG even though you asked it to save as a geotiff. Have you customized your environment/configs at all? What happens if you pass filename="{name}_{start_time:%Y%m%d_%H%M%S}.tif"?

FYI I edited your comment to put your traceback in three backticks (`) on each end. It shows up better for code literals than quoted blocks (>).

[dtlloyd]

Thanks for amending my comment and helping me out with this problem.

I don't recall customising the environment or configs at all. I am running SatPy in a Jupyter notebook. Could that be a problem?

As for passing a filename="{name}_{start_time:%Y%m%d_%H%M%S}.tif" to the save_datasets method, I get the same error as before. If I change the datatype argument to dtype=np.int32, I get no error, but the output .tif image is blank (uniform pixel value).

Perhaps if I told you what I am trying to do, you may be able to suggest a work around. I have resampled the Sentinel 3 dataset to a new map projection. I now want to either save one of the Sentinel 3 bands as a geotiff and work on that, or write the same Sentinel 3 band to a numpy array and work on that. Is there a different way I can do that in SatPy?

[djhoese]

Oh, so you just want to work with the data? You can use the xarray DataArray object directly with scn['band_name']. In the long run it would be best to learn how to use xarray, but if you need to get a numpy array:

my_data_arr = scn['band_name']
my_np_arr = my_data_arr.compute().data

The DataArray objects in satpy have dask arrays underneath. Calling compute computes these dask arrays which makes them numpy arrays. Using .data gets that numpy array from the DataArray. Note that you are now holding on to all of that data in memory which can fill up quickly.

I'd like to still investigate this issue you're having but I'm not exactly sure where to go from here. I'll try some stuff locally.

[dtlloyd]

Fantastic advice @djhoese. That's exactly what I wanted to do. I am not sure how I missed that in the documentation.

If there is anything more I can do to help you debug this issue please let me know.

[djhoese]

@dtlloyd Do you have rasterio installed? Can you do import rasterio? What version of trollimage and satpy are you using?

Edit: I'm not able to reproduce this with ABI data.

[dtlloyd]

satpy version 0.13.0
trollimage version 1.7.0

I have rasterio installed but not imported. It was imported in another notebook open at the same time as the satpy one I am working from.

I am afraid that I can't replicate the bug now. I have tried restarting the kernel a couple of times and I no longer get an error when I run newscn.save_datasets(writer='geotiff', dtype=np.float64, enhance=False).
I honestly don't know what's happened there. Sorry for the bother.

[djhoese]

I'm calling this closed. If someone has this same issue again then please open a new issue and reference this. Thanks.