Make temerature scene

[gersmit]

Hi @djhoese,

I now have nice images in color, beautifull, grid, coast, borders, deco, points, everything oke.
But small question about things i try.
In matplotlib.pyplot i can make a temprature images of Europe, i try this with an composite routine but canot find how to make it.
Can you explane to me how to get same result with a composite -> local_scn.save_dataset

THIS IS THE TEST FROM MAKE TEMPERATURE MATPLOTLIB - PYPLOT _ CARTOPY:
#brightness - temp IR_039
print('Creeer warmte kaart Europa')
import cartopy.crs as ccrs
ir = 'IR_039'
coor_box(-10,35,25,58)
local_scn = scn_cropped
crs = local_scn[ir].attrs['area'].to_cartopy_crs() #ccrs.PlateCarree()
fig, ax = plt.subplots(figsize=(11,7)) # start plot
ax = plt.axes(projection=crs) #ccrs.LambertConformal()'wereld'
my_datak = local_scn[ir]
my_datac = local_scn[ir]-273.15 # plot data cel
my_datac.plot.imshow(transform=crs)
ax.coastlines()
ax.gridlines()
ax.axis('tight') # maak border small
fig.tight_layout() # alleen border grafiekbar
fig.savefig(data_img_warmte + 'Europa_temperature ' + str(file_name) + '.png')
#plt.show()

Thanks in advance for an answer
@ger

[djhoese]

So the above code works or it doesn't work? You said you tried to make it work with an RGB composite? The above code seems to be using IR_039 which is not a composite.

Additionally, what errors are you getting? It is difficult to tell what is going wrong without any error messages or without telling me exactly how it isn't working or what is unexpected.

[gersmit]

The code works fine for me but what you say, with the IR_039.
What i want to no is how to make the same output but then from the compositor.
How can i use or read the data from the IR_039 in a composite.

The errors are different because i make misstakes.

When i read in de docs the working of the composite it is logical.
I try it and use it now.
But to start to amke a composite as i ask, the logiacl is not there.
Maybe you give me a hint how to start to get the IR_039 and use the data for the color.

[djhoese]

Do you maybe need to use the get_enhanced_image function shown in this part of the tutorial: https://github.com/pytroll/tutorial-satpy-half-day/blob/master/notebooks/05_composites.ipynb

It would help to see the code you are using to generate the composite so I can know exactly what code could be used to plot it. If you have the composite DataArray object you should be able to do:

data_arr = local_scn[my_composite_name]
from satpy.writers import get_enhanced_image
img = get_enhanced_image(data_arr)
enhanced_data_arr = img.data
enhanced_data_arr.plot.imshow(rgb='bands', transform=crs)

[gersmit]

@djhoese
To make things clearer i make AreaDefinition.

In example "Remote Sensing II" (area_id = "Griechenland) they use llx = 10E5 as coordinate.
I will use coordinate in meter : llx = -895620.3449964523, but then the image is match to small north to south.
Or i wil use degrees.

Can you explane to me what is the value 10E5 and how can i calculate these to degrees.
And why is the image when i nuse meters to small (N/S) and with the same values in yaml oke.

[djhoese]

In example "Remote Sensing II"

What is this? Is that part of this tutorial?

Can you explane to me what is the value 10E5 and how can i calculate these to degrees.
And why is the image when i nuse meters to small (N/S) and with the same values in yaml oke.

I can't make a guess without seeing the actual code and YAML you are using. 10E5 is the same as 1000000:

$ ipython
Python 3.7.6 | packaged by conda-forge | (default, Mar 23 2020, 23:03:20)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.19.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: 10e5
Out[1]: 1000000.0

For future questions, please include at least the code you used and any error messages it produced. If something is unexpected, tell me why. If something doesn't work, tell me in what way it doesn't work. If an image doesn't look right, show me a screenshot of the image.

[gersmit]

In the attachement is a temperture image, with the mouse i can find the value by coordinate/pixel (ru).
Is there a way to get the coordinate/pixel value by Satpy script so i can get the value to a file.

[djhoese]

If you know the lon/lat coordinate or x/y coordinate or the col/row of the pixel, then yes.

row = 5
col = 5
scn[composite].values[row, col]

If you know the lon/lat then you can do:

lon = 1.25
lat = 5.33
area = scn[composite].attrs['area']
col, row = area.get_xy_from_lonlat(lon, lat)
# then do the same as above

[gersmit]

Thanks, 👍👍

[gersmit]

@djhoese, is there for getting RGB color value the same.
Btw, is there something where i can find examples to get parameter values fron scenes?

[djhoese]

RGB should work the same exact way. Try it and let me know if it doesn't work.

Btw, is there something where i can find examples to get parameter values fron scenes?

I'm not sure what you mean. Isn't that what I showed you above?

[gersmit]

yes, the examples you write to me

[gersmit]

Because i wil get the pixel color to see if the agreement is the same with the temperature, I try :
local_scn[ir].RGB[col, row]
Error =
Traceback (most recent call last):
File "/home/pi/a-satpy-testen/test_overlay_maak_sneeuw_kaart.py", line 509, in
local_scn[ir].RGB[col, row]
File "/home/pi/.local/lib/python3.7/site-packages/xarray/core/common.py", line 229, in getattr
"{!r} object has no attribute {!r}".format(type(self).name, name)
AttributeError: 'DataArray' object has no attribute 'RGB'

And a question, i place attachement 2 days ago, the white circle upper left come from creating the area.
Is there a possibilty to remove these circles?

[djhoese]

AttributeError: 'DataArray' object has no attribute 'RGB'

Where did you get the code to do local_scn[ir].RGB[col, row]? Is that in a tutorial somewhere?

You should be accessing your composite directly local_scn[your_rgb_name].values[row, col]. NOTE: You had row, col backwards as well.

Is there a possibilty to remove these circles?

Try passing a larger radius_of_influence to your Scene.resample call. Something like radius_of_influence=50000.

[gersmit]

Dear @djhoese,
Getting better every time.
Question:
I have AreaDefinition with the image i wil use.
area_id = "Europa"
description = "Europa Mercator-Projektion"
proj_id = "Europa"
proj_dict = {"proj": "merc", "lat_ts": 38, 'lon_0': 25}

width = 1000 # width of the result domain in pixels
height = 675 # height of the result domain in pixels

llx = -32E5 # projection x coordinate of lower left corner of lower left pixel
lly = 27E5 # projection y coordinate of lower left corner of lower left pixel
urx = 10E5 # projection x coordinate of upper right corner of upper right pixel
ury = 68E5 # projection y coordinate of upper right corner of upper right pixel

area_extent = (llx,lly,urx,ury)
area_def = AreaDefinition(area_id, proj_id, description, proj_dict, width, height, area_extent)

But how can i make it so that i can use lat/lon coordinates insteat of llx = - 32E5.
When it is possible it is easyer to use with calculations.
How can i calculate the distance of one pixel in the plot from Matplotlib.

Thanks for an answer and for being patient to inform us

[djhoese]

How can i calculate the distance of one pixel in the plot from Matplotlib.

I'm not sure about matplotlib, but if you have an AreaDefinition you can do area_def.pixel_size_x or area_def.pixel_size_y to get the pixel size in the projection space.

For using degrees instead of meters, you can look at using the create_area_def function: https://pyresample.readthedocs.io/en/latest/geometry_utils.html#areadefinition-creation. Something like this might work:

area_def = create_area_def(area_id, proj_string, width=1000, height=675, area_extent=(lon_min, lat_min, lon_max, lat_max), units='degrees')

NOTE: That the "merc" (mercator) projection is still in meters even though you defined your area in degrees. That's just the way projections (and the PROJ C library) work.

[gersmit]

Works well, thanks.
I have now made several AreaDefinition, how do I determine the height and width in pixels so that everything remains in format.
Is there a formula for that?

[djhoese]

in format

What do you mean in format? Please look at the pyresample documentation that I linked above. There is also a page specifically for AreaDefinition objects and how to create them. In my above code I specify width and height so I'm not sure what else you need.

[gersmit]

I think if the area changes, the height and width also change

[djhoese]

It depends how you create the area. With the create_area_def, I think (as long as you don't change the width and height provided) if you change the area_extent it won't effect the width and height.

[gersmit]

I create a scene as shown in the appendix earlier.
I am using width 1000pix and height 675 pix.
Now I create AreaDefinition lon_min = -12 lat_min = 28 lon_max = 28 lat_max = 64.
With the tool col, row = area.get_xy_from_lonlat (lon, lat) everything works fine and the place names made with plt.text are in the right place.
Now I change the AreaDefinition to a smaller size (Netherlands) and the place names are no longer correct. Height and width is the same.
Can you explain to me how this can be solved, where to adjust. I can't find what the line area = scn [composite] .attrs ['area'] means. What is meant by .attrs ['area']. I can not fond the variable area

[djhoese]

Code! I need to see code. You say:

Now I change the AreaDefinition to a smaller size (Netherlands) and the place names are no longer correct. Height and width is the same.

How are you doing this? If you print this area print(area_def) what does the output look like?

Can you explain to me how this can be solved, where to adjust. I can't find what the line area = scn [composite] .attrs ['area'] means. What is meant by .attrs ['area']. I can not fond the variable area

The Scene object (scn) is a container for xarray DataArray objects. So scn[composite] is a DataArray object. DataArray objects have metadata defined in a .attrs dictionary. So that's scn[composite].attrs is a dictionary of metadata about the composite. Satpy fills in this information when you load the composite. One of those pieces of metadata is the 'area'. This area is the AreaDefinition defining the geolocation of the data. So scn[composite].attrs['area'] is the AreaDefinition for composite. Does that make sense?

Edit: If you get an error when doing scn[composite].attrs[area], what is it?

[gersmit]

@djhoese, i am sorry to tel that it is oke now.
I forgot to change a value in a new routine.
my apologies but you explaine to me another issue with AreaDefinition.

[gersmit]

Dear @djhoese,
For working with HRV from native i use now virtual Raspberry pi Buster.
Lot of space and when i use the composite

from satpy.composites import RealisticColors
compositor = RealisticColors("realcols", lim_low=85., lim_high=95.)
composite = compositor([local_scene['VIS006'],
... local_scene['VIS008'],
... local_scene['HRV']])
It is running.
It take long time to build an image that i save to a map.
But withoud eny errorr there is not a file saved.
Can you explane to me what i do wrong.
Other writer perhaps.
I use:
ocal_scn.save_dataset('sat_image', meteo_file ,
overlay = {'coast_dir': '/home/pi/sat_data/shape_data_root/',
'overlays': {'grid': grid,
'coasts': coast, 'borders': borders,
'rivers': rivers, 'lakes': lakes, 'points':points}},
decorate = {'decorate': deco}, fill_value=0)

This is working good with other composite.

[djhoese]

Could you please provide all the code? HRV is bigger than the other datasets. A Raspberry Pi is not powerful. It will take time and memory.

[gersmit]

I use now RPI in Vmware running in Windows.
Very powerfull and enoufh memeory.
The above routine i use simply with:

from satpy.writers import to_image
img = to_image(composite)
img.invert([False, False, True])
img.stretch("linear")
img.gamma(1.7)
img.show()
With other composite no problem.
And what i say, everything runs, i takes his time to transleate the HRV but when it is allright it wil not save as i save normaly.
There are no errors.
Mabey you have an test routine that use HRV and show and write to image.

[djhoese]

Could you please provide all the code?

ALL the code

I can not guess at what could be going wrong if I don't know everything you are doing (imports, scene creation, call to load, call to resample, everything). If it is too much code, then make a smaller version that still produces the error.