Surface Pressure Overlays with PyTROLL/Satpy


Ernst Lobsiger
 

Dear All,

here is some more fine stuff I intend to make accessible in a feature update of the PyTROLL/SatPy Starter Kit. Projections is where Satpy shines due to the underlying Proj libraries.
The 5th Corona wave is hitting Switzerland. Around Christmas I predict lockdown with plenty of time at home. This here is from the NOAA Ocean Prediction Center. It's today 12:00 UTC.

Cheers,
Ernst


Christian Peters
 

Ernst,

great stuff. Looking forward to the new PyTROLL/SatPy Starter Kit update. 
Overlays are really a great enhancement for received satellite data and amateur weather forecaster. 
Stay healthy.

Regards,

Christian 

Am 03.12.2021 um 20:44 schrieb Ernst Lobsiger via groups.io <ernst.lobsiger@...>:

Dear All,

here is some more fine stuff I intend to make accessible in a feature update of the PyTROLL/SatPy Starter Kit. Projections is where Satpy shines due to the underlying Proj libraries.
The 5th Corona wave is hitting Switzerland. Around Christmas I predict lockdown with plenty of time at home. This here is from the NOAA Ocean Prediction Center. It's today 12:00 UTC.

Cheers,
Ernst <Meteosat-11-20211203-SLO-1200-airmass-opcae-wc.jpg><GOES16-20211203-SLO-1200-airmass-opcaw-wc.jpg><GOES17-20211203-SLO-1200-airmass-opcpe-wc.jpg><Himawari-8-20211203-SLO-1200-airmass-opcpw-wc.jpg>


Andreas Mueller
 

Ernst,

can only echo that. Looking forward to the new Pytroll version as well.

Andreas


Graham Woolf
 

Hi Ernst

Me too  - it would be great if you could also show how you arrive at the pixel size for the overlay for the geographical area chosen

I hope the maths isnt too complicated

Regards

Graham


Ernst Lobsiger
 

On Sat, Dec 4, 2021 at 01:41 AM, Graham Woolf wrote:
Me too  - it would be great if you could also show how you arrive at the pixel size for the overlay for the geographical area chosen
Graham,

I already posted point by point what has to be done. The pixel size is mainly given by the original chart you use und how you make it ready with IM as an overlay.
The true problem can be to figure out the map projection used. As far as the "Bracknell Chart" is concerned I finally got a useful answer from UK MetOffice today:

<crs-definition name="EPSG:9999" type="proj4" proj4-string="+proj=stere +lat_0=90 +lon_0=35w +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"/>

Needless to say that this is exactly what I figured out with a ruler, a triangle and a pocket calculater after I got a dummy automatic MetOffice response two weeks ago..

 
Regards,
Ernst





Ernst Lobsiger
 

On Sat, Dec 4, 2021 at 01:41 AM, Graham Woolf wrote:
I hope the maths isnt too complicated
Graham and All,

Here is a little competition for those of you who do not just want to click and pray the Windows way:

I attached a forcast chart easy.gif for 2021/12/05/00:00. You have plenty of time to get that working.
You can overlay that in a couple of easy steps on any GEOS16 composite. Show us your results tomorrow.

1) You should be capable to find out that easy.gif is purly BW (Paint 3D) and has a size of 1728 x 1100.

2) You must have ImageMagick convert.exe installed. ** DON'T TRY C:\Windows\System32\convert.exe! **.
   Users of my PyTROLL/Satpy Starter Kit 3.0 already have it under C:\EMCtools\exefiles\convert.exe

3) Convert easy.gif to an overlay 32Bit *.png making all white transparent using convert.exe from 2)
   Open a CMD window and change directory where you have your easy.gif (maybe even desktop) then type:
   C:\EMCtools\exefiles\convert easy.gif -transparent white PNG32:easy-bw.png
   If you rather want it white on transparent you have to negate the input chart first. Just try this:
   C:\EMCtools\exefiles\convert easy.gif -negate -transparent black PNG32:easy-wb.png

4) Now you have tho funny looking overlays that you can use if Satpy makes you a 1728 x 1100 image
   that has the same projection as your original chart. It's easily seen that the central Meridian
   is about 65W and the projection most probably is a Cylindrical Mercator (I was too lazy to test).
   So make an entry in areas.yaml : Name it easy, type any descriptive text, try projection merc,
   use Ellipse WGS84, give it a central meridian like this  lon_0: -65.0  , you can also define a
   parameter   lat_0: 20.0   which is about the center of your image. I no lat_0 is given it defaults
   to lat_0 = 0.0. Now under shape give it the pixel size of easy.gif. Finally some "try and error":
   You have to indicate the area extent in meters. It is measured from the point (lon_0, lat_0).

   https://proj.org/operations/projections/index.html               (for users that want to know more)

5) Now make any kind of GOES16 image with area='easy' and compare (especially the grid) with easy.gif.
   If everything fits you use a final statement to apply your overlay. For this, study some posts starting here:

   https://groups.io/g/MSG-1/message/32242                           (read this and my following answers)

   You can further improve the fit comparing the two grids and by fine tuning the 'easy' area_extent.


Good luck,
Ernst (the teacher)



David J Taylor GM8ARV 🏴󠁧󠁢󠁳󠁣󠁴󠁿 🇪🇺
 

On 04/12/2021 14:44, Ernst Lobsiger via groups.io wrote:
Graham and All,

Here is a little competition for those of you who do not just want to click and
pray the Windows way:
The OS on which the software is run should not matter, and it's a confusion to
refer to the OS in this way. You will find much Linux software which no longer
requires compilation and arcane commands to make it work. What I do see,
though, and have experienced myself, is the infinite upgrade loop in Linux.
People reporting that they have needed to upgrade a program for some software
to work, only to find that some other software no longer works and needs to be
upgraded, forcing yet more updates. I see this less in Windows.

Many people are not interested in the internals of a software set, nor in
complex programs which do everything but are not easy to use. They do just
want to "click and go", and get the weather images and data they need. Whether
that's on a Windows system, an Android or Apple tablet, or on a Raspberry Pi
400 with a good display doesn't matter. For other folk, the internals may be
of as much interest as the data itself!

Cheers,
David
--
SatSignal Software - Quality software for you
Web: https://www.satsignal.eu
Email: david-taylor@blueyonder.co.uk
Twitter: @gm8arv


Ernst Lobsiger
 

David,

the headline of this thread is "Surface Pressure Overlays with PyTROLL/Satpy".
It's for people that use my Kit and thus have shown that they do want and can
do more than "click and go". It's true that by the crazy idea to bring GNU/Linux
to the desktop it has been windonized more and more. There is a growing base of
GNU/Linux users that want "click and go" because they do not and cannot do more.

My receipt is for the few interested ones and works on most if not any OS. The
difference is e.g. that GNU/Linux users that have IM installed just type convert
while under Windows there is a tiny naming conflict with a Windows convert.exe.

This is also for people that know that they cannot afford (or even buy at all)
a "click and go" program that has the flexibility and power of PyTROLL/Satpy.
And PyTROLL/Satpy will be ready and still *free* on day one we get MTG1 data.
My Kit has been run under GNU/Linux, Windows 10, Mac OS and on Raspberry Pi.


Cheers,
Ernst

P.S. It's the same thing with DVB-S2 receivers: There is no "click and go" to make a
TBS-5927 box work well or a TBS-6909X work at all under Windows. But both are perfect
receivers under GNU/Linux. If you want Windows "click and go" you buy a Novra router.


Ernst Lobsiger
 

On Sat, Dec 4, 2021 at 06:44 AM, Ernst Lobsiger wrote:
   you can also define a  parameter   lat_0: 20.0   which is about the center of your image. I no lat_0 is given it defaults  to lat_0 = 0.0.
Sorry, this is my bad. The Mercator projection has no parameter lat_0 or does set it 0 whatever you set it in area 'easy' . The area_extent is always measured from the equator (lon_0, lat_0=0).

https://proj.org/operations/projections/merc.html

There is a parameter  lat_ts  which means your map will have a true scale (1m on earth = 1m on map) at this latitude. If you omit it it defaults to 0 (true scale along the equator which is easiest).

Cheers,
Ernst




Graham Woolf
 

Hi Ernst

Here is my image

As you can see the difference between the lines of latitude at 10 and 30 degrees are greater than the difference between the 20 degree lat . How do I get the 10 and 30 deg lines on the overlay close to the grid on the image without moving the 20 deg lat line on the overlay

Regards

Graham


David J Taylor GM8ARV 🏴󠁧󠁢󠁳󠁣󠁴󠁿 🇪🇺
 

On 06/12/2021 09:32, Graham Woolf wrote:
Hi Ernst

Here is my image

As you can see the difference between the lines of latitude at 10 and 30
degrees are greater than the difference between the 20 degree lat . How do I
get the 10 and 30 deg lines on the overlay close to the grid on the image
without moving the 20 deg lat line on the overlay

Regards

Graham
Graham,

When setting the projections for GeoSatSignal, I would change the latitude span
for an offset like you show. This might also require changing the parallel.

Cheers,
David
--
SatSignal Software - Quality software for you
Web: https://www.satsignal.eu
Email: david-taylor@blueyonder.co.uk
Twitter: @gm8arv


Ernst Lobsiger
 

On Mon, Dec 6, 2021 at 01:32 AM, Graham Woolf wrote:
Here is my image

As you can see the difference between the lines of latitude at 10 and 30 degrees are greater than the difference between the 20 degree lat . How do I get the 10 and 30 deg lines on the overlay close to the grid on the image without moving the 20 deg lat line on the overlay
Graham,

well done! You need just some fine tuning. Things will be easier if you make the grid drawn by satpy red to easy distinguish it from the overlay grid.
Here is a general receipt. The process is iterative so you must make a couple of images until everything fits. Finally you can turn the satpy grid off.

1 If a red meridian at the left side of the image is left of the corresponding
  white meridian of the overlay we must go left (and vice versa).
  We must shift the x value of lower_left_xy  to the left.

2 If a red meridian at the right side of the image is right of corresponding
  white meridian of the overlay we must go right (and vice versa).
  We must shift the x value of upper_right_xy  to the right.

3 If a red parallel near the top of the image is above the corresponding
  white parallel of the overlay we must go up (and vice versa).
  We must shift the y value of upper_right_xy  further up.

4 If a red parallel near the bottom of the image is below the corresponding
  white parallel of the overlay we must go down (and vice versa).
  We must shift the y value of lower_left_xy  further down.


Applying this simple receipt you must go up on top and down at the bottom. The center parallel will not move much. And you have to go left on the left side and also (a little bit less) left on the right side.

Good luck,
Ernst


Graham Woolf
 

Hi Ernst

Thanks

I have got the lines a bit closer but not sure where to go now. They seem OK in the middle but not at the edges
 
this is the x-y I am using

  area_extent:
    lower_left_xy: [-4489999.0,      0.0]
    upper_right_xy: [3888888.0, 3000000.0]

Regards

Graham


Ernst Lobsiger
 

On Mon, Dec 6, 2021 at 05:20 AM, Graham Woolf wrote:
I have got the lines a bit closer but not sure where to go now. They seem OK in the middle but not at the edges
 
this is the x-y I am using
Graham,

you have to make finer steps. This process must converge. The steps you make must be finer and finer down to maybe 2000m at the end !

Look at the red 90W meridian. This one was left of the white 90W and following the receipt you went left.
But now the red 90W is right of the white 90W you must go back right (you went too far left in the first step, go back right half).

The same applies to the red 30N parallel. It was slightly above the white 30N and you were supposed to go slightly up.
But it seems you went down because the red 30N is now way up and even out of sight. It's now the red 20N that tells you to go up.

It's easy to mix up x and y, left and right, plus and minus. You have to concentrate and double ckeck.
Maybe you make a # commented version of the area_extent and note there your improvements.
Then if you make the next adjustment you see where you come from.

Good luck,
Ernst


Ernst Lobsiger
 

Graham and All,

here is the solution for this litte exercise:

easy:
  description: NOAA radiofax forcast chart showing wind and wave height in the western atlantic
  projection:
    proj: merc
    ellps: WGS84
    lon_0: -65.0
  shape:
    height: 1100
    width:  1728
  area_extent:
    lower_left_xy: [-4230000.0, -655000.0]
    upper_right_xy: [3790000.0, 4420000.0]

Lessons learnt: The Satpy parser is picky. Make 2 spaces indents and at least 1 space after the colon ":"  .
The safest way is to copy an existing entry in areas.yaml, rename it 'easy' and change the parameters.

The test overlay was choosen for its simplicity (pure black and white, no antialiasing) and taken from here:

https://www.nhc.noaa.gov/tafb_latest/atl24_latestBW.gif

Maybe you find more usesful stuff to overlay on this site. Attached a heavily reworked overlay from OPC.
It shows storm BARRA (HARRY) from yesterday. Some sites in the UK might have to repoint their dishes.

Cheers,
Ernst


Christian Peters
 

Ernst and all,

another Overlay made with Pytroll/Satpy and with nice overlay code form Ernst: 
It’s the IR_108 inverted channel, masked with NWC-CT data on top of blue marble background with an DWD overlay made with a script of Ernst, glued with IM. 
Maybe available in the next Starter Kit from Ernst…!? ;-)

Regards,

Christian
 

Am 08.12.2021 um 10:56 schrieb Ernst Lobsiger via groups.io <ernst.lobsiger@...>:

Graham and All,

here is the solution for this litte exercise:

easy:
  description: NOAA radiofax forcast chart showing wind and wave height in the western atlantic
  projection:
    proj: merc
    ellps: WGS84
    lon_0: -65.0
  shape:
    height: 1100
    width:  1728
  area_extent:
    lower_left_xy: [-4230000.0, -655000.0]
    upper_right_xy: [3790000.0, 4420000.0]

Lessons learnt: The Satpy parser is picky. Make 2 spaces indents and at least 1 space after the colon ":"  .
The safest way is to copy an existing entry in areas.yaml, rename it 'easy' and change the parameters.

The test overlay was choosen for its simplicity (pure black and white, no antialiasing) and taken from here:

https://www.nhc.noaa.gov/tafb_latest/atl24_latestBW.gif

Maybe you find more usesful stuff to overlay on this site. Attached a heavily reworked overlay from OPC.
It shows storm BARRA (HARRY) from yesterday. Some sites in the UK might have to repoint their dishes.

Cheers,
Ernst

<Meteosat-11-20211207-SLO-1200-overview-opcae-wc.jpg>


Alan Curnow
 

That looks most impressive Christian. Thanks to Ernst too for taking the time developing all these extra features for Satpy/Pytroll. 

Barra was quite blowy Ernst, but I tightened the clamps on my big dish a bit more after it moved a few weeks ago, and it's currently still pointing in the right direction, according to the S/N :)
Alan 


Ernst Lobsiger
 

Dear All,

in the latest GEO Newsletter No 72 on page 21 there is an article
"Overlay MetOffice Surface Pressure Charts over Meteor M2 Images".

http://leshamilton.co.uk/GEO/PDF/geoq72.pdf

But you can hardly fit a Polar Stereographic onto an UTM, which is
basically a Transverse Mercator. As both projections are conformal,
(retaining angles and small shapes) you maybe can, by rotating and
scaling, fit the UK more or less. Then Spain and Iceland are way off.

There is only two ways to do this mathematically right:

1) You can, if your system and software allows, reproject the
   satellite data the same way MetOffice made the MSLP chart.
2) If you cannot reproject your satellite image, then you can
   reproject the MSLP chart as overlay of your satellite image.

Maybe you have never seen 2) but this can be done in any GIS that
deserves this name (e.g. Q-GIS is free and works in almost any OS).

Last sunday Christian and myself solved method 2) in PyTROLL/Satpy.
This means we can now reproject the Bracknell chart or any of the DWD
charts to any other area (except for dynamic areas like the 'omerc_bb').
This does take less than 1 minute and it can easily be automated.

I attach two absolutely extreme reprojection examples. I fit the
Bracknell chart and a DWD Type n chart (Stereographic and Conic)
to a NOAA Ocean Prediction Center (OPC) Mercator of Atlantic East.
Think of these charts as printed on a thin rubber sheet. You then
distort this sheet to exactly fit to your target projection lat/lon grid.
Reprojecting raster data comes at a price of lower quality though.

Cheers,
Ernst


Ernst Lobsiger
 

Dear users of MeteorGIS,

I have no Quadrifilar Helix, no SDR stick and no MeteorGIS here. But I doubt that MeteorGIS can easily do what I demonstrated below. And if you reproject a raster map in Q-GIS you must exactly know the UTM version MeteorGIS uses. If you really want a Bracknell chart overlay on your Meteor images here is my two cents how it can possibly be done right (no guarantee). AFAIK MeteorGIS accepts shape files to reproject. Here is a link that hints how you can use some GIS functionality of MeteorGIS (unfortunately the connection is rather unstable):

https://usradioguy.com/layering-fire-and-hotspot-date-on-meteor-sat-imagery

I have almost no experience with Q-GIS. But I made a quick V 3.16 install and did this:

1) I imported one of my Bracknell overlays that I exported to GEO tiff with PyTROLL/Satpy.
    If you only have a *.gif or *.png you will have to specify the projection parameters.
    Here is what I got from the UK MetOffice:     https://groups.io/g/MSG-1/message/32274

2) Q-GIS has GDAL under the hood that allows you to polygonize (vectorize) a raster map.

3) This vectorized layer can easily be reprojected to EPSG_4326 (lon +/- 180°, lat +/- 90°).


You can then export this layer as shape file(s) and use it in MeteorGIS just like coast lines.
If someone is interested in trying this I can send him the shapes (5MB) or publish them here.

Cheers,
Ernst


I attach screen shots of the vectorized Bracknell chart (note the outlined fonts!).
and of the final shape file in Plate Carrée "projection" (just lon/lat in degrees).


Ernst Lobsiger
 

Dear users of MeteorGIS,

I contacted Les Hamilton and he sent me an offline version of a MeteorGIS. It took some time and a couple of e-mails to get that running. I can only make one single image and do not have a useful version of the Bracknell chart for that date. So the clouds and the pressure lines do not match. It took more time to find out that MeteorGIS only undestands single point and single polyline shapes. The documentation of MeteorGIS is very sparce. But we did it as explained in the last post.

Cheers,
Ernst