#!/usr/bin/python3 # -*- coding: utf-8 -*- #************************************************************************* # # NOAA-20 JPSS (Joint Polar Satellite System) sensor viirs M-Bands # https://directory.eoportal.org/web/eoportal/satellite-missions/n/noaa-20 # https://www.wmo-sat.info/oscar/satellites/view/noaa_20 # # Orbit: Sun-synchronous app. LTAN=13:30 (Local Time on Ascending Node) # # ***** This script takes an argument YYYYmmDD ***** # # It only takes daylight passes (raw swath images are upside down) #------------------------------------------------------------------------- # # CH-3123 Belp, 2021/01/17 License GPL3 (c) Ernst Lobsiger # #************************************************************************* # # File Pattern : SVMC_{platform_shortname}_d{start_time:%Y%m%d_t%H%M%S%f}_e{end_time:%H%M%S%f}_b{orbit:5d}_c{creation_time:%Y%m%d%H%M%S%f}_eum_ops.h5 # Example name : SVMC_j01_d20190411_t0224358_e0226003_b07218_c20190411023355000121_eum_ops.h5 # Hugo's ruler : 0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012 # copied EMCV : 0 1 2 3 4 5 6 7 8 9 10 # EUMETCast : These files arrive as above in channel NPP-2 (Basic Service) # ANNOTATION : This script uses ImageMagick (IM) convert for annotation # # I need from datetime import datetime, timedelta from glob import glob import math import sys import os # From PyTroll/Satpy we take from pyorbital.orbital import Orbital from satpy import Scene, MultiScene # Needed for grid labels that have their own color setting (see above grid definition) # This did segfault for a while and has then been fixed by Dave (may or may not work)! # import aggdraw # Why to hell is it not working? # from satpy.utils import debug_on # debug_on() # Minimal command line parameter test usage = sys.argv[0] + ' expects YYYYmmDD as CLI parameter:\n' \ + 'See comment concerning EUMETCast in the script header.' if len(sys.argv) != 2: sys.exit(usage) if len(sys.argv[1]) != 8: sys.exit(usage) # Date YYYYmmDD Dat=sys.argv[1] # Get substrings Yea = Dat[:4] Mon = Dat[4:6] Day = Dat[6:] # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv # ********** ADAPT THESE PARAMETERS TO YOUR PERSONAL NEEDS *********** # Edit 1-3 parameter(s) below according to your file layout: # My files are in a /MountPoint/Channel/YYYY/mm/DD structure # No trailing / because this will make timestamps unreadable segdir = '/srv/rec_0/NPP-2' + '/' + Yea + '/' + Mon + '/' + Day # Windows toodrp = 'C:', datdrp = 'D:' / GNU/Linux use homedir toodrp = '/home/eumetcast' datdrp = '/home/eumetcast' # Available composites: See files ../satpy/etc/composites/visir.yaml and viirs.yaml # ********************************************************************************* # ['ash', 'cloudtop_daytime', 'dust', 'false_color', 'fire_temperature', # 'fire_temperature_39refl', 'fire_temperature_awips', 'fire_temperature_eumetsat', # 'fog', 'ir108_3d', 'ir_cloud_day', 'natural_color', 'natural_color_sun_lowres', # 'natural_with_night_fog', 'night_fog', 'ocean_color', 'overview', 'snow_age', # 'snow_lowres', 'true_color', 'true_color_lowres', 'true_color_lowres_crefl', # 'true_color_lowres_land', 'true_color_lowres_marine_tropical', 'true_color_raw'] composite = 'natural_with_night_fog' # Please invent 'comp' as short form of too long 'composite' used for IM annotation # Sorry, this kind of abbreviation cannot be automatically set in a meaningful way! comp = 'natural_wnf' # Center coords (° North and ° East) and range ° of interest # Westminster lat = 51.228 lon = -0.532 ran = 30.0 # Isle of Man # lat = 54.228 # lon = -4.532 # ran = 20.0 # Switzerland # lat = 46.88 # lon = 8.15 # ran = 20.0 # Area definition must consider above coords and ran for best results # Do not use 'swath' or SatPys dynamic areas 'omerc_bb' and 'laea_bb' # Choose your area # area = 'eurol' # area = 'isleofman' area ='westminster' # area = 'switzerland' # Configure your individual overlays, either True/False or 1/0 do work ADDgrid = True ADDcoasts = True ADDborders = True ADDrivers = True ADDpoints = True # OVRCache speeds up static overlays. These are stored in your .../EMCdata/cache # Never use this with homebrew or SatPy's dynamic areas like omec_bb and laea_bb # With caching you will sometimes have to delete *.png OVRcache files manually ! OVRcache = False # ******** TOUCH THE CODE BELOW ONLY IF YOU KNOW WHAT YOU DO ********* # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # If you set tlefil = '' it will be downloaded from Internet! tlefil = toodrp + '/EMCtools/pppconfig/my_TLE_file.txt' tmpdir = datdrp + '/EMCdata/tmpdirs/xnoaa20' imgdir = datdrp + '/EMCdata/images/NOAA-20' # Dir must exist try: os.chdir(tmpdir) except: sys.exit('Cannot change to tmp directory ...') # NOAA-20 is an afternoon satellite heading north # LTAN = 13:30 (sometimes also called ECT asc) # LTDN = 01:30 (sometimes also called ECT desc) d1 = int(13. + 30./60. - lon/15. - 5.5) d2 = int(13. + 30./60. - lon/15. + 6.5) # Degrees to radian dr = math.pi/180.0 # Use Spherical Trigonometry: # In Memoriam Prof. Paul Wild sinlat = math.sin(lat*dr) coslat = math.cos(lat*dr) cosran = math.cos(ran*dr) # Initialize goodfiles = [] goodtimes = [] scenes = [] cosmax = 0 offset = len(segdir)+1 orb = Orbital('NOAA-20', tlefil) # Glob by date, assume that interesting segments are all after midnight files = glob(segdir + '/' + 'SVMC_j01_d' + Yea + Mon + Day + '*.h5') # Find satellite position from file name for fname in files: # These are times at segment begin Hou = fname[offset+20:offset+22] Min = fname[offset+22:offset+24] Sec = fname[offset+24:offset+26] # Day crossings between d1 and d2 if int(Hou) > d1 and int(Hou) < d2: # Segments are 85 - 86 seconds long, let's advance to the center dt = datetime(int(Yea), int(Mon), int(Day), int(Hou), int(Min), int(Sec)) + timedelta(seconds = 43) # Get longitide, latitude, altitude of sat (slon, slat, salt) = orb.get_lonlatalt(dt) # Use Nautical Triangle for arc distance cosarc = sinlat * math.sin(slat*dr) + \ coslat * math.cos(slat*dr) * math.cos((lon-slon)*dr) # Within search range, maybe more than 1 orbit if cosarc > cosran: goodfiles.append(fname) goodtimes.append(int(dt.timestamp())) # We only want the best orbit if cosarc > cosmax: cosmax = cosarc timmax = int(dt.timestamp()) if goodfiles == []: sys.exit('Sorry, no good files found ...') # Get useful info back from timmax dt = datetime.fromtimestamp(timmax) orbmax = orb.get_orbit_number(dt) telmax = dt.strftime('%H:%M UTC') tfname = dt.strftime('%H%M') # The blend function actually stacks the passes one after the other numpass = 0 # Shingle passes from the east towards the center timfirst = min(goodtimes) while timmax - timfirst > 1500: passfiles = [] n = len(goodtimes) - 1 while n >= 0: if goodtimes[n] - timfirst < 1500: passfiles.append(goodfiles[n]) goodfiles.pop(n) goodtimes.pop(n) n -= 1 timfirst = min(goodtimes) numpass += 1 passfiles = sorted(passfiles) print('Pass ', numpass) # DEBUG: for n in range(0, len(passfiles)): print ('%2d --> ' % (n+1), passfiles[n]) scenes.append(Scene(filenames = passfiles, reader = 'viirs_compact')) # Shingle passes from the west towards the center timlast = max(goodtimes) while timlast - timmax > 1500: passfiles = [] n = len(goodtimes) - 1 while n >= 0: if timlast - goodtimes[n] < 1500: passfiles.append(goodfiles[n]) goodfiles.pop(n) goodtimes.pop(n) n -= 1 timlast = max(goodtimes) numpass += 1 print('Pass ', numpass) # DEBUG: passfiles = sorted(passfiles) for n in range(0, len(passfiles)): print ('%2d --> ' % (n+1), passfiles[n]) scenes.append(Scene(filenames = passfiles, reader = 'viirs_compact')) # Now there should just be the center pass left. We stack it on top for best possible image quality. # This might not always work as expected as cosarc is only calculated for the center of the segment. goodfiles = sorted(goodfiles) numpass += 1 print('Pass ', numpass) # DEBUG: for n in range(0, len(goodfiles)): print ('%2d --> ' % (n+1), goodfiles[n]) scenes.append(Scene(filenames = goodfiles, reader = 'viirs_compact')) print('Found', numpass, 'NOAA-20 pass(es) in range to stack (and possibly image depending on chosen area)') print('The closest pass (with max elevation at POI) is Orbit number', orbmax,' peaking at', telmax) # Does also work with just 1 pass (if ran has been set too small!) # For MultiScene (Experimental) see the latest satpy documentation mscn = MultiScene(scenes) mscn.load([composite]) new_mscn = mscn.resample(area, radius_of_influence = 5000) #, reduce_data = False) stacked_scn = new_mscn.blend() # Next is line by Christian Peters, thanks a lot! # try:except: added after Graham Woolf's problems try: # This is for real (multichannel) composites layers, height, width = stacked_scn[composite].shape except: # In case of a single channel 'composite' height, width = stacked_scn[composite].shape # Autoscale font_size (fs) of grid labels fs = int((height + width) / 100) # my_font = aggdraw.Font('red', datdrp + '/EMCdata/fonts/DejaVuSerif.ttf', # opacity = 255, size = fs) # aggdraw.Font takes only 4 arguments # Font color can either be specified 'white' or (R, G, B) = (255, 255, 255) # .. 'font': my_font} below may or may not work anymore now or in the future. # If a direct font definition does not work, use the inline definition below: # New interface to pycoast + POIs have been integrated 2020, thanks Martin !! my_grid = {'major_lonlat': (10,10), 'minor_lonlat': (2, 2), # Opacity 0 will hide the line, values 0 ... 255 EL 'outline': (255, 255, 255) , 'outline_opacity': 255, 'minor_outline': (200, 200, 200),'minor_outline_opacity': 127, 'width': 1.5, 'minor_width': 1.0, 'minor_is_tick': False, 'write_text': True, 'lat_placement': 'lr', 'lon_placement': 'b', 'font': datdrp + '/EMCdata/fonts/DejaVuSerif.ttf', 'font_size': fs} # minor_is_tick: False draws a line not ticks, True does ticks # label placement l,r,lr for latitude and t,b,tb for longitude # Colors are specified as strings (140 HTML color names) or as triple (R, G, B) 8Bit. # Best 'resolution' is automatically choosen according to area_def, no need to touch. # Coast also has keys ... 'fill': 'green', 'fill_opacity': 50} for land mass filling. # 'level' defaults to 1 but should be choosen at least 2 for rivers (1 = 0 in GSHHG). my_coasts = {'outline': (255, 255, 0), 'width': 1.5, 'level': 1} # , 'resolution': 'i'} # level 1 .. 4 my_borders = {'outline': (255, 0, 0), 'width': 1.0, 'level': 1} # , 'resolution': 'h'} # level 1 .. 3 my_rivers = {'outline': ( 0, 0, 255), 'width': 1.0, 'level': 3} # , 'resolution': 'f'} # level 1 .. 11 # List of POIs (Station, Cities, etc. ...) my_poi_list = [((13.4052, 52.5214), 'Berlin'), (( 2.3522, 48.8566), 'Paris'), (( 0.1278, 51.5074), 'London'), ((-3.1880, 55.9539), 'Edinbourgh'), ((-6.2598, 53.3514), 'Dublin'), (( 7.4989, 46.8918), 'Belp')] # You may find other True Type Fonts (*.ttf) in /usr/share/fonts/truetype/... # Windows 10 PRO users may find more True Type Fonts in C:\Windows\Fonts\ ... my_points = {'font': datdrp + '/EMCdata/fonts/DejaVuSerif.ttf', 'font_size': 20, 'points_list': my_poi_list, 'symbol': 'circle', 'ptsize': 20, 'outline': 'black', 'fill': 'red', 'width': 5.0, 'outline_opacity': 255,'fill_opacity': 255, 'box_outline': 'navy', 'box_linewidth': 3.0, 'box_fill': 'gold', 'box_opacity': 255} # Cache is useful for speedup. Does not work with dynamic 'omerc_bb' or 'laea_bb' !! my_cache = {'file': datdrp + '/EMCdata/cache/DAY/' + composite, 'regenerate': False} # Setup indivdual overlays dictionary my_overlays = {} if ADDgrid: # Add key grid my_overlays['grid'] = my_grid if ADDcoasts: # Add key coasts my_overlays['coasts'] = my_coasts if ADDborders: # Add key borders my_overlays['borders'] = my_borders if ADDrivers: # Add key rivers my_overlays['rivers'] = my_rivers if ADDpoints: # Add key points my_overlays['points'] = my_points if OVRcache: # Add key cache my_overlays['cache'] = my_cache # Enhance image with overlays and save *.png to temporary directory # Note the slightly different method to save (e.g. save_datasets!!) stacked_scn.save_datasets(filename = 'NOAA20.png'.format(area = area), overlay = {'coast_dir': datdrp + '/EMCdata/gshhg-shp', 'overlays' : my_overlays}) # 2^6=32 combinations should do # **************************************************************** # NEW POST PROCESSING: WORKS WITH ALL SCRIPTS AND PROJECTION AREAS # **************************************************************** # ImageMagick (IM) postprocessing, added directory logos in EMCdata # logo1 = datdrp + '/EMCdata/logos/GNU-Linux-Logo-Penguin-SVG.png' logo1 = datdrp + '/EMCdata/logos/NOAA_1200x1200.png' logo2 = datdrp + '/EMCdata/logos/PyTROLL_400x400.jpg' # Use IM to annotate (maybe also -equalize) and convert to *.jpg # www.imagemagick.org/Usage/crop/#splice or .../annotating/#labeling # This is the syntax used -draw "text OffsetX,OffsetY 'Hello World!'" # This annotates to the left, numbers depend on original image size # Colors can be {dark,light}blue, etc. or rgb(R,G,B) with \( escaped def magick_string(EffHeight): """ It seems ImageMagick under Windows cannot cope with multiline text that includes \n newlines. We have to draw every line separately. This scales all pixel values for IM according to the image height. Should make it easier to use annotation and logos at the left side. This has been optimized using MSG4 full scan with 3712x3712 pixels. """ OptHeight = 3712 OptSplice = 640 OptLogo1s = 540 OptLogo1x = 50 # 40 for 'GNU-Linux-Logo-Penguin-SVG.png' OptLogo1y = 40 # 60 for 'GNU-Linux-Logo-Penguin-SVG.png' OptLogo2s = 530 OptLogo2x = 55 OptLogo2y = 55 OptPoints = 90 fac = EffHeight/OptHeight EffSplice = str(int(fac*OptSplice+0.5)) EffLogo1s = str(int(fac*OptLogo1s+0.5)) EffLogo1x = str(int(fac*OptLogo1x+0.5)) EffLogo1y = str(int(fac*OptLogo1y+0.5)) EffLogo2s = str(int(fac*OptLogo2s+0.5)) EffLogo2x = str(int(fac*OptLogo2x+0.5)) EffLogo2y = str(int(fac*OptLogo2y+0.5)) EffPoints = str(int(fac*OptPoints+0.5)) legend = [' \'DATE: (DAY)\'\" ', ' \''+Yea+'/'+Mon+'/'+Day+'\'\" ', ' \'TIME:\'\" ' , ' \''+telmax+'\'\" ', ' \'SOURCE:\'\" ' , ' \'EUMETCast\'\" ', ' \'SERVICE:\'\" ' , ' \'Basic\'\" ', ' \'CHANNEL:\'\" ' , ' \'NPP-2\'\" ', ' \'RECEIVER:\'\" ' , ' \'TBS-6909X\'\" ', ' \'SATELLITE:\'\" ' , ' \'NOAA-20\'\" ', ' \'SENSOR:\'\" ' , ' \'viirs\'\" ', ' \'COMPOSITE:\'\" ' , ' \''+comp+'\'\" ', ' \'AREA:\'\" ' , ' \''+area+'\'\" '] multidraw = '' yoff = fac*470 for n in range(0, len(legend)): if (n % 2) == 0: xoff = fac*40 yoff += fac*150 else: xoff = fac*60 yoff += fac*100 multidraw += '-draw \"text +'+str(int(xoff+0.5))+'+'+str(int(yoff+0.5))+legend[n] cmdstr = 'convert NOAA20.png '\ + '-gravity west -background rgb\(245,245,245\) -splice '+EffSplice+'x0 '\ + logo1+' -gravity northwest -geometry '+EffLogo1s+'x'\ + EffLogo1s+'+'+EffLogo1x+'+'+EffLogo1y+' -composite '\ + logo2+' -gravity southwest -geometry '+EffLogo2s+'x'\ + EffLogo2s+'+'+EffLogo2x+'+'+EffLogo2y+' -composite '\ + '-pointsize '+EffPoints+' -gravity northwest '+multidraw\ + imgdir+'/'+'NOAA-20-stack-'+Dat+'-DAY-'+tfname+'-'+composite+'-'+area+'.jpg' return(cmdstr) # *Image Magick* os.system(magick_string(height)) # IM will still be working when this python script returns (subprocess.call?)