Script to generate Overshooting Top

Post Reply
User avatar
Posts: 1
Joined: Fri Jun 05, 2020 9:23 am

Script to generate Overshooting Top

Post by bonyseptian »

Any one have idea to make a script in McIDAS V for generating the Overshooting Cloud Top from IR 10.8 (Himawari8) netcdf dan NWP Tropopause Temperature netcdf data following Bedka (2010) algorithm:

1) Read 10.7 IR Window (IRW) channel BT and NWP Tropopause Temperature

2) Identify pixels with IRW BT < 215 K and NWP tropopause temp and filter in a list

3) Starting from the coldest BT in the list, ensure that no cold pixels are within 15 km of each other (each IR pixel have resolution 2x2 km). Pixels that satisfy these criteria are called candidate overshooting pixels

4) Sample the surrounding anvil at an 8 km radius in 16 directions around the candidate overshooting pixel. 8 km is chosen to sample the region outside an OT, since OTs are generally < 15 km in diameter. We don't want to sample an OT when computing mean surrounding anvil BT

5) Pixels at least 6.5 K colder than surrounding anvil are considered as overshooting center pixels

6)Within a search box centered on the OT center pixel, find pixels at least 50% colder than the surrounding anvil mean BT. It is used to identify all "TOP" and not just one pixel on "Top"
User avatar
Posts: 938
Joined: Mon Nov 15, 2010 5:57 pm

Re: Script to generate Overshooting Top

Post by bobc »


This is the same information that I sent you through email earlier. I'll post it here as well so it is available to others.

In short, I believe it is possible for McIDAS-V to do everything you are looking for, but there is no straight-forward/easy way of doing it. This would require scripting. If you haven't done so already, I'd encourage you to visit the McIDAS-V Documentation page where you'll find two scripting tutorials (Basic and Advanced Scripting). These will give an introduction to basic scripting functionality in McIDAS-V followed by some more in-depth analysis of data objects.

I put together a script that you could possibly consider as a starting point. Continue reading if you want to give it a try:

I'm attaching a script (ot.txt - you can rename it to that you can try running through the Jython Shell in McIDAS-V. Note that the script begins with a comment of a probe() function you'll need to add to your Jython Library. Also, you'll likely need to change the date/time passed through loadADDEImage and the URL passed through loadGrid because the servers/catalogs are realtime and don't contain data that is that old, especially the satellite data. Essentially, this script does the following:

  1. Loads a GOES-16 band 13 temperature image (small subsetted domain). You can change this up to another server/dataset you have access to if you are interested in a different satellite.
  2. Loads GFS tropopause temperature data at the same timestep as the satellite data. This is a GFS dataset over the United States. You might want to consider a different global/regional model over your area of interest.
  3. Performs a mask on both the satellite and gridded data so only pixels with temperatures less than 215K are included in the data objects. You can individually display these if you uncomment the createLayer lines.
  4. Creates a data object containing only pixels where both the satellite and grid mask thresholds are satisfied.
  5. Uses various functions/methods described in the Advanced Scripting tutorial to extract lat/lon values from the masked data object created in step 4.
  6. Uses the probe() function (from the Advanced Scripting tutorial and the attached script) to get the temperature value at each lat/lon in the masked data object.
  7. Creates a string of "temperature, lat, lon" at each pixel, and removes all of the "missing" values since these are from pixels that don't satisfy the mask criteria.

At this point, you may be able to work with various javadocs for VisAD and McIDAS-V to explore things a bit further. Please note that this is just the way that I began attempting the process described in your email, it's possible that there are different ways of doing this. Perhaps others on the forums have done this type of work in McIDAS-V before and would be able to give you some more pointers.

Bob Carp
McIDAS User Services

Here is the content of ot.txt:

Code: Select all

def probe(field, lat, lon):
    from import GridUtil
    from visad.georef import EarthLocationTuple
    altitude = 0.0
    loc = EarthLocationTuple(lat, lon, altitude)
    valueAtLocation = GridUtil.sampleToReal(field, loc, None)
    return valueAtLocation

# load satellite imagery data source
addeParms = dict(
    server = '',
    dataset = 'RTGOESR',
    descriptor = 'CONUSC13',
    unit = 'TEMP',
    #size = 'ALL',
    coordinateSystem = LATLON,
    size = (500, 500),
    location = (22.8, -90),
    day = '2020-06-03',
    time = '06:01:13',
satData = loadADDEImage(**addeParms)

# load gridded data source
gridParms = dict(
    filename = 'dods://',
    field = 'Temperature_tropopause',
    time = '2020-06-03 06:00:00Z',
gridData = loadGrid(**gridParms)

# create masked data source of satellite imagery data
# and display only where temperatures are <= 215K
maskDataS = mask(satData, '=<', 215, 1) * satData
#maskLayerS = activeDisplay().createLayer('Image Display', maskDataS)
#maskLayerS.setLayerLabel('Satellite Temps =< 215K')

# create masked data soruce of gridded data
# and display only where temperatures are <= 215K
maskDataG = mask(gridData, '=<', 215, 1) * gridData
#maskLayerG = activeDisplay().createLayer('Image Display', maskDataG)
#maskLayerG.setLayerLabel('Grid Tropopause Temps =< 215K')

# display masked satellite data only at pixels where both
# satellite and grid mask conditions are satisfied
maskDataSG = mask(satData, '=<', 215, 1) * maskDataG
maskLayerSG = activeDisplay().createLayer('Image Display', maskDataSG)
maskLayerSG.setLayerLabel('Satellite and Grid Tropopause Temps =< 215K')

# get lat lons
domainSet = getDomainSet(maskDataSG)
latLons = getLatLons(domainSet)
lats = latLons[0]
lons = latLons[1]

# use probe function for each lat/lon point
dataList = []

numPts = len(lats)
for x in range(0, numPts):
    dataList.append(probe(maskDataSG, lats[x], lons[x]))

# create list of dataVal, lat, lon
newList = []
for x in range(0, numPts):
    newList.append('%s - %s - %s' % (dataList[x], lats[x], lons[x]))
# modify newList to not include any points that have dataVal of missing
newList2 = [x for x in newList if not x.startswith('missing')]
Rename to
(2.4 KiB) Downloaded 75 times
Post Reply