Getting brightness temperature values according to different colors in MSG-SEVIRI ash RGB composite

[Thaize Baroni]

Hello,

I'm new to satpy and I wonder if somebody could help me. There are 2 things that I need to do using the High Rate SEVIRI Level 1.5 Image Data - MSG - 0 degree.

1) I'm using SEVIRI's ash composite to distinguish volcanic ash from ice/water clouds. Here is the code I wrote to visualize the ash composite:

from satpy import Scene
import matplotlib.pyplot as plt
from satpy.writers import get_enhanced_image
from pyresample import geometry
from pyresample import create_area_def
import numpy as np

dateien = ["/home/data/seviri/W_XX-EUMETSAT-Darmstadt,VIS+IR+HRV+IMAGERY,MSG4+SEVIRI_C_EUMG_20181224130009.nc"]

files = {"seviri_l1b_nc" : dateien}
scn = Scene(filenames=files)
scn.load(['ash'])

area_def = create_area_def('my_area', {'proj': 'lcc', 'lon_0': 16, 'lat_0': 37.5, 'lat_1': 47.5},

                 width=3000, height=3000,
                 area_extent=[12, 36, 18, 39], units='degrees')
local_scn = scn.resample(area_def)

crs = local_scn["ash"].attrs["area"].to_cartopy_crs()

fig = plt.subplots(figsize=(10,10))
ax = plt.axes(projection=crs)
ax.coastlines(resolution="10m", color="white")
gl=ax.gridlines(xlocs=range(12,18,1),ylocs=range(36,39,1), draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False

img = get_enhanced_image(local_scn['ash'])
img_data = img.data
img_data.plot.imshow(rgb='bands', vmin=0, vmax=1)

plt.title("MSG at " + local_scn.attrs["start_time"].strftime("%Y-%m-%d %H:%M"))

My question is: It is possible to "access" the different colours from the ash composite? If yes, how?

For example, in the ash composite, the ash-contaminated pixels have a reddish colour. It is possible to get all the brightness temperature values only from the pixels with reddish colours? If yes, how?

2) I need to estimate the minimum brightness temperature at 10.8 in an area of 9x9 pixels centered over a specific point (latitude 37.7510 and longitude 14.9934).

How would I do that? 

Regards,

Thaize

[David Hoese]

Hi Thaize,

1. Your `img_data` variable is the RGB color of the ash RGB. If you have
a way to detect what pixels you consider to have enough red. If you can
get that then you can create a mask or index array to use on any other
channel you load. For example:

scn.load(['ash', 'IR_108'])
... your existing code ...

img_data = img.data
# True where red pixels are greater than 0.8 and green and blue are less
than 0.2
red_mask = (img_data[:, :, 0] > 0.8) & (img_data[:, :, 1] < 0.2) &
(img_data[:, :, 2] < 0.2)
# .data is a dask array, we do .compute to get a numpy array
bt_data = local_scn['IR_108'].data.compute()
bt_for_red_pixels = bt_data[red_mask]

This is a simple approach and there are many optimizations that could be
made. Additionally, it may be easier or more accurate for you to figure
out what makes the ash RGB red-ish and do those calculations yourself to
come up with a better red_mask than what I've done below (use the bands
that make up the ash RGB).

2. The first thing to try might be to use the Scene's crop method:
https://satpy.readthedocs.io/en/stable/api/satpy.html#satpy.scene.Scene.crop

If that doesn't do what you want you could look at the methods for the
AreaDefinition (scn['IR_108'].attrs['area']) in pyresample for how to
convert lon/lat to the col/row in your array:
https://pyresample.readthedocs.io/en/latest/api/pyresample.html#pyresample.geometry.AreaDefinition.lonlat2colrow

That would give you the nearest index (row, column) to the lon/lat you
have which you could then use to slice the array:

your_square = scn['IR_108'][my_row - 3:my_row + 4, my_col - 3: my_col + 4]

I think that would work.

You could also look at making your own AreaDefinition. See
https://pyresample.readthedocs.io/en/latest/geometry_utils.html. You can
then resample to that AreaDefinition and the result would be the region
you want.

Hope that helps.

Dave

> --
> You received this message because you are subscribed to the Google
> Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pytroll+u...@googlegroups.com
> <mailto:pytroll+u...@googlegroups.com>.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/pytroll/b739d91a-92d9-45c2-a2f9-1ff9b70630d2n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/b739d91a-92d9-45c2-a2f9-1ff9b70630d2n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Thaize Baroni]

Hi Dave,

Thank you for your help!

I like the idea of creating a mask for the red pixels. However, when I tried to do what you suggested, I got:

bt_for_red_pixels = bt_data[red_mask]
IndexError: too many indices for array: array is 2-dimensional, but 3 were indexed

So, I masked the colours individually, doing the following:

img = get_enhanced_image(local_scn['ash'])
img_data = img.data

####green
img_data_green=img_data[1, :, :]
img_data_green=np.asarray(img_data_green)
img_data_green[img_data_green >0.4]=np.nan ###masking all the non ash pixels
green=pd.DataFrame.from_records(img_data_green)
######red
img_data_red=img_data[0, :, :]
img_data_red=np.asarray(img_data_red)
img_data_red[img_data_red <0.8]=np.nan ###masking all the non ash pixels
red=pd.DataFrame.from_records(img_data_red)

red.update(green)
img_data_final=np.asarray(red) 
bt108min=np.asarray(local_scn['IR_108']) 
df_bt108min=pd.DataFrame.from_records(bt108min)
df_bt108min=df_bt108min[red.notnull()] ### keeping only ash-related pixels

I know it is longer and dumbier, but it worked!

Now, let's say I want to use yellow-ish pixels as well. How would I do that?

Regards

[David Hoese]

Ah sorry, I must have given you the wrong index locations. It should
have worked the way I gave it to you, but you got it to work so no problem.

For yellow, it sounds like you have all the building blocks you need.
Determine what "yellow" means to you and for your use case and do the
masking accordingly.

Dave

> https://groups.google.com/d/msgid/pytroll/e3b39602-3cd0-4b61-99b9-45b9436ac312n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/e3b39602-3cd0-4b61-99b9-45b9436ac312n%40googlegroups.com?utm_medium=email&utm_source=footer>.