Issue when including redband in Rayleigh.get_reflectance

[aschae11]

Problem description

When I run my code with red band included I get the traceback below, however, if I call the function excluding the red band it runs as expected. I noticed at the bottom of the 'rayleigh.py' script the test function doesn't include the red band.

I have checked that all of my arrays going into the function are of the same shape, so I am rather confused why it doesn't like the shape of the red band.

Actual Result, Traceback if applicable

Traceback (most recent call last):
File "/Volumes/GoogleDrive/My Drive/GOES_PROJECT/merge_full_disk_for_color_with_gamma.py", line 107, in
blue_cor = ray.get_reflectance(sun_alt, satz, ssadiff, 'C01', red)
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pyspectral/rayleigh.py", line 228, in get_reflectance
redband = from_array(redband, chunks=redband.shape)
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dask/array/core.py", line 2327, in from_array
token = tokenize(x, chunks)
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dask/base.py", line 606, in tokenize
return md5(str(tuple(map(normalize_token, args))).encode()).hexdigest()
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dask/utils.py", line 415, in call
return meth(arg)
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dask/base.py", line 745, in normalize_array
data = hash_buffer_hex(x.copy().ravel(order='K').view('i1'))
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/ma/core.py", line 3046, in view
output = ndarray.view(self, dtype)
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/ma/core.py", line 3271, in setattr
self._mask.shape = self.shape
ValueError: total size of new array must be unchanged

Versions of Python, package at hand and relevant dependencies

Python 2.7.15, Pyspectral 0.8.2

Thank you for reporting an issue !

[adybbroe]

Hi @aschae11
Can you post the code producing the above traceback?
As well as version of numpy and dask.

[aschae11]

Current numpy is 1.11.0 and dask is 0.18.2.

Below is the condensed version of my script that still produces the issue.

from future import division
import numpy as np
from netCDF4 import Dataset as nc
from pyspectral.rayleigh import Rayleigh
from pyorbital.astronomy import get_alt_az, sun_zenith_angle
from pyorbital.orbital import get_observer_look
import datetime

date = '20180828'
time = '1600'

data_dir = '/Users/schaefer/Google Drive File Stream/My Drive/GOES_PROJECT/GOES/FULL/'
iname_blue = "GOES16_EAST_FULLDISK_BAND_01_"+date+"T"+time+".nc"
iname_red = "GOES16_EAST_FULLDISK_BAND_02_"+date+"T"+time+".nc"

ifile = nc(data_dir+iname_blue)
dims = ifile.dimensions
nx = dims['x'].size
ny = dims['y'].size

blue = ifile.variables['Sectorized_CMI'][:,:]
sat_lat = ifile.satellite_latitude
sat_lon = ifile.satellite_longitude
sat_alt = ifile.satellite_altitude
ftime = ifile.observation_time
ftime = datetime.datetime.strptime(ftime, '%Y%m%dT%H%M')
ifile.close()

ifile = nc(data_dir+iname_red)
red = ifile.variables['Sectorized_CMI'][:,:]
ifile.close()

lon = np.arange(-165+180/1809, 15, 180/1809)
lat = np.arange(-90+180/1809,90, 180/1809)

lats, lons = np.meshgrid(lat, lon)

sun_alt, sun_ang = get_alt_az(ftime,lons,lats)
sun_ang = sun_ang % 360
sata, satel = get_observer_look(sat_lat, sat_lon, sat_alt/1000, ftime, lons, lats, 0)
sata = sata % 360
satz = 90 - satel
ray = Rayleigh('GOES-16','abi')
ssadiff = np.absolute(sun_ang - sata)
ssadiff = np.minimum(ssadiff, 360-ssadiff)

blue = np.minimum(blue,1.0)
blue = np.maximum(blue,0.0)
red = np.minimum(red, 1.0)
red = np.maximum(red, 0.0)

blue_cor = ray.get_reflectance(sun_alt, satz, ssadiff, 'C01', red)

[adybbroe]

Thanks @aschae11
I am not familiar with these files, but are you sure the red and blue data arrays have the same shape?
The GOES red and blue bands does not have the same resolution.
I have files like
OR_ABI-L1b-RadF-M3C01_G16_s20172631730408_e20172631741175_c20172631741218.nc
and
OR_ABI-L1b-RadF-M3C02_G16_s20172631730408_e20172631741175_c20172631741209.nc

I can try test your script if you post the data somewhere I can fetch them.

May I ask why you do not use SatPy to do all this. SatPy reads ABI data and does atm correction with Pyspectral and resample data using Pyresample.

Did you see the examples like this one for Himawari:
https://github.com/pytroll/pytroll-examples/blob/master/satpy/ahi_true_color_pyspectral.ipynb
?

You can do the same with ABI, just replace the reader name with abi_l1b

[aschae11]

@adybbroe,

You probably don't recognize the name because they are files I created from data coming in from our LDM feed.

Since I am working with the full disk images, the red and blue bands do have the same resolution (6 km), though you are right that they are different, 0.5 km and 1 km respectively, when you get to the CONUS and MESO files.

Google drive links to the two files:
https://drive.google.com/file/d/1CCJ2Tlvpk_9jKHqEyVuAa5zISsot2yIZ/view?usp=sharing
https://drive.google.com/file/d/19zvOfgmn02Cqc42Ge6h_hbDn78ckr5B4/view?usp=sharing

I am doing all of this "free hand" because I am honing my skills and building my portfolio for reference when I start applying for jobs. I may try satpy to see how well it runs, but I am curious if it will handle my files as expected. Most of the algorithms that have been used (sim_green) were dug out of satpy by me or others so I am sure it is a good package, though I am also curious the processing speed. Some of the code seemed bulky as I was digging through it.

Thanks for looking into this.

[djhoese]

though I am also curious the processing speed. Some of the code seemed bulky as I was digging through it.

A couple things about this: It looks like you are using SCMI files for ABI data. SatPy does not have the ability to read SCMI files directly right now but it is something that I want to (and will need to) add in the future. Any help would be much appreciated.

About satpy's performance and/or bulky-ness: code bulky-ness does not necessarily say anything about the performance of a piece of software. The SatPy code has also evolved over many years (starting from the old mpop package). SatPy also handles a lot of different cases and lets users do a lot of different things so it needs some complexity to the internals to handle these cases. If you have suggestions for things that didn't make sense to you or are very confusing, let me know or make an issue on the satpy repository asking for better documentation for a particular component. We have a lot of examples so I would hope that the highest level interfaces can handle most things.

As for performance, satpy depends on xarray and dask which allow us to have very fast processing code. The true color ABI image below was made with satpy in 8 minutes on my machine. This image uses the rayleigh correction from pyspectral, has been sharpened using a ratio sharpening technique, has been corrected for sun angle, and has been reviewed by multiple scientists to produce the best looking results we know how to make. Again 8 minutes with a maximum of ~12GB of memory!

Edit: Additionally the benefits of working together (i.e. satpy) greatly out weigh any benefits of working alone.

[aschae11]

Hi @djhoese,

I am working with the Sectorized CMI, It seemed goofy to convert to radiance to just end up converting that to RGBs, so I stayed with SCMI. It does make trying to remove the Rayleigh scattering more challenging.

I agree that bulkiness of code doesn't say anything about performance, it just makes it a little more challenging to completely follow along coming at it from an outside point of view. I completely understand the need for complexity since you are trying to make the package applicable to many situations.

What ratio sharpening technique did you use for the example image? I assume both the ratio sharpening and sun-angle corrections are embedded within satpy too?

Currently, I am trying to turn images over a little faster than 8 minutes. I will eventually be creating full disk and CONUS images, so I need turnover to be quicker. Currently, my image processes in under one minute, on both my laptop and school machine (8 and 9 GB RAM, respectively). This compute time includes the rayleigh correction (to an extent, because I cheated and converted the output radiances from pyspectral into SCMI values) and applying a simple gamma filter to brighten things up. The output example is attached below. What file did you use to create your image? Looks like a different position than GOES16.

Edit: Is it even worth running the ratio sharpening on the full disk since all of the bands are the same resolution? It makes sense for CONUS and MESO sectors with the red band resolution difference but is it making a difference at the even 6km resolution of full disk?

[pnuu]

@djhoese is using the full resolution data with 1 km resolution. That is, 36 times more data compared to your 6 km version. Could be that he made the example image using GOES-17 preliminary test data.

The sun zenith angle correction uses PyOrbital to calculate the illumination angles, see sstpy.composites.SunZenithCorrectorBase class.

[djhoese]

@aschae11 Regarding your first paragraph regarding SCMI and radiances, what do you mean converting to radiances? It is true that the ABI L1b files do start as radiances but are calibrated to BTs and reflectances. The BT and reflectances in the SCMI files that I'm familiar with are produced with the same calibration calculations. The abi_l1b reader in satpy can produce Refl/BT or radiance data, but the reflectance data is what is used for making RGBs so there is no need to "unconvert" data back to radiances.

RE bulky code: Understood. Just keep in mind that if we can find a way to work together then we all benefit.

As @pnuu mentioned, but I'll correct, my image is made with the full disk full resolution data from GOES-16 ABI. So that means the red band is 500m and the lower resolution bands required are 1km resolution but replicated to a 500m resolution. I guess you could consider this 144x more data in the final result. The data I used for the above image is from late 2017 when GOES-16 was still in the "TEST" position (center longitude -89.5) which is why it looks a little "off" geographically.

The ratio sharpening is something that Liam Gumley (SSEC/CIMSS, where I work) suggested I try where the red band (500m) has every 4 pixel block averaged then replicated to the full resolution. This "averaged" band is then compared with the value of the full resolution red then multiply that ratio by the other lower resolution bands. In most cases it seems to improve the overall look of the image. Since we are dealing with full resolution (500m and 1km) reflectance bands in this RGB the sharpening does make sense. You are right that it makes less sense if the red band has been decimated to a lower resolution from the start.

I have one question about your method, what do you mean when you say that you convert pyspectral radiances in to SCMI values? Pyspectral (at least the Rayleigh correction part) can be passed reflectances directly which is what I've seen typically stored in the tiled SCMI files that are provided to AWIPS. Instead of a gamma correction we've been playing with a couple log-like scalings including this one and a multi-point linear enhancement (very similar to log/sqrt) here.

I'll see if I can find some lower resolution data somewhere or manually decimate it and see what kind of performance I can get.

Edit: About satpy's code complexity, it is definitely something we want to work on so if you have suggestions we'd love github issues or pull requests. If you have any non-SCMI data you could also try some of the existing examples on ABI L1B data and see what type of performance you get for full resolution data.

[aschae11]

@djhoese

Regarding your first paragraph regarding SCMI and radiances, what do you mean converting to radiances?

Using the CMI data directly to create the 3 bands works great, however, when I start applying the different corrections, the CMI data becomes more difficult to work with.

When I would try and use the CMI data with pyspectrals' rayleigh scattering, I would end up with output values much greater than the CMI data. My interpretation of the pyspectral function was that including the red band allows it to be scaled to the data being used, however the original issue relates directly back to why it probably was not working as expected. I did find the information needed (buried in the CMIP documentation) of how to go from CMI to the Radiances output from GRB (i.e. the 'L1b' files). Once converted to the Radiances, the rayleigh scattering function works as expected, even when including the red band. Similarly, working with the L1b data, the function works as expected (it seems to really darken the images but that my just be because of my lack of knowledge for the rest of the image processing.

I think it would be fairly easy to implement the conversion to be able to use the CMI data, though I think working with the L1b files makes life much easier due to the added projection data included in the files.

I have one question about your method, what do you mean when you say that you convert pyspectral radiances in to SCMI values?

I noticed the SCMI values are usually between 0 and ~1.05, while the Radiance values from the L1b files range ~0-800ish? If pyspectral is expecting the values of the latter then subtracting 24 makes sense, however subtracting 24 from the former nullifies the image.

I'll see if I can find some lower resolution data somewhere or manually decimate it and see what kind of performance I can get.

The Thredds test server has a lot of the 6km decimated files (under goes 16) , along with the full size data (under GRB16).

Instead of a gamma correction we've been playing with a couple log-like scalings including this one and a multi-point linear enhancement (very similar to log/sqrt) here.

I like the first scaling however I don't understand the second one, It will take deeper digging to understand it.

I would like to pick your brain about the data flow within Satpy and what format (range or values) it is at different stages.

[djhoese]

@adybbroe Can correct me, but the rayleigh scattering in pyspectral should have worked regardless. That said, it may be dependent on the reflectances provided being in ~0-120% range rather than the typical albedo values of ~0-1.2. This may be why your radiance input/output looks more like what you expect, but if pyspectral is expecting reflectances then your output is probably still incorrect.

FYI I've just added a SCMI reader to satpy but I haven't officially merged it to master yet. You can see the status of this feature here: pytroll/satpy#416
It should work for any SCMI files that act as one giant image (no tiles). I did it in the effort to combine satpy with metpy which you can see here: https://github.com/pytroll/pytroll-examples/blob/master/satpy/SatPy-MetPy-Siphon%20ABI%20L1B%20over%20THREDDS.ipynb

I noticed the SCMI values are usually between 0 and ~1.05, while the Radiance values from the L1b files range ~0-800ish? If pyspectral is expecting the values of the latter then subtracting 24 makes sense, however subtracting 24 from the former nullifies the image.

I'm having trouble following everything you are saying. I think you're forgetting to give me some details. As mentioned above reflectances are usually either in albedo (units: 1) between 0 and ~1.20 or as a percentage (units: %) between 0 and ~120. You say "24" but where did you get that number from? Regardless I think it is safe to assume that pyspectral expects reflectances in percentages but may accept the albedo units too (not sure).

The Thredds test server has a lot of the 6km decimated files (under goes 16) , along with the full size data (under GRB16).

Again, I think you forgot to mention some details about your process. This is the first time you've brought up a THREDDS server.

Please join our slack team if you haven't already so we can chat more about the data flow of satpy. There is a link for inviting yourself to the team under the "getting in touch" section on the pytroll website: http://pytroll.github.io/

[adybbroe]

@aschae11
I tested your code above under Python 3.6.4 It ran without problems here.
However, I have found one issue. The sun angles returned by get_alt_az from pyorbital are in radians, not degrees!

Also, the atmospheric correction in PySpectral works with reflectances (in percent), not radiances.
See https://pyspectral.readthedocs.io/en/master/rayleigh_correction.html

Hope this helps?
I have kind of lost track of what all your issues are, if several?

You are welcome asking on Slack, I think it will be easier trying to understand how we can help you best that way.

-Adam