How do I customise Default Background Maps?

Post any questions, ideas, or topics related to Jython and Python scripting.
User avatar
bobc
Posts: 987
Joined: Mon Nov 15, 2010 5:57 pm

Re: How do I customise Default Background Maps?

Post by bobc »

Hello,

Thanks for the explanation of why operating on a layer before displaying it would be ideal for you. I wrote up inquiry 2991 to evaluate adding this functionality, though I don't have an estimate on when this would be done.

As you noticed, there is no ability to load data through a script that emulates the ability of the capabilities of the HYDRA chooser. I wrote up evaluating this as inquiry 2992.

One difference between the images is that through scripting you are loading the BRIT (brightness, 0-255 scale) calibration, and the HYDRA chooser uses the TEMP (temperature) calibration for emissive bands. I would suggest changing:

Code: Select all

unit='BRIT'
to

Code: Select all

unit='TEMP'
in your script and see if that improves things. I'm a bit unsure about the white stripes in your data loaded through loadADDEImage. I have a few granules I'm testing with and I'm not seeing that problem. Could you please post your MOD02 file to our anonymous ftp server for me to take a look at? Here are instructions of how to do that through a command prompt/terminal window:

Code: Select all

ftp ftp.ssec.wisc.edu
user: anonymous
password: email
bin
prom
cd pub/incoming
put filename
bye
Please let me know the name of the file once you post it and I'll give it a look.

Thanks,
Bob
User avatar
m.yunus.a.m.
Posts: 18
Joined: Wed Aug 18, 2021 3:55 am

Re: How do I customise Default Background Maps?

Post by m.yunus.a.m. »

Hi Bob,

I have uploaded the file, it's name is MOD021KM.A2016142.1020.061.2017326015248.hdf. I downloaded it from https://ladsweb.modaps.eosdis.nasa.gov/ ... /2016/142/. I'm using Windows 10 on i5-4460 if it helps. No GPU, my GPU just bricked recently.

Regarding using the unit='TEMP', it works great. After manually adjusting the range to match the range loaded through HYDRA, the image is similar to the one loaded by HYDRA, except for the white bands. It jumps around when the Texture Quality is adjusted. Several Texture Quality options would greatly distorts the image though.

That being said, how does HYDRA determines the range to be used in the Image Display, and how do I replicate it in scripting?

Yunus
Attachments
adde-temp.png
User avatar
bobc
Posts: 987
Joined: Mon Nov 15, 2010 5:57 pm

Re: How do I customise Default Background Maps?

Post by bobc »

Hello,

When loaded through the HYDRA chooser, McIDAS-V sets the min/max enhancement values based on the min/max of the data in the scene. When loaded through ADDE, McIDAS-V uses a parameter default to display the data with a range of 180 - 320K.

You can see statistical information about the data object returned from loadADDEImage using McIDAS-V's describe() function. For example, in the Jython Shell you can run:

Code: Select all

print(describe(Image1))
You'll see that the min is listed as 0.0 and the max is 486.29. It's common for the minimum value to list as 0.0. Normally, the max value returned from describe() would actually match the max enhancement value used in a display from the HYDRA chooser. It turns out that there are actually 40 data points with this 486.29 value, which I verified with a Grid Table display. There looks to be a bug using the South Pole projection in that you see the missing gaps of data right along the dateline. I believe we already have an open inquiry to handle some similar issues with the dateline.

One workaround you can do to mask out these extraneously high/low values would be to use McIDAS-V's maskWithinRange() function to set values outside of an expected range to be missing. For example:

Code: Select all

maskedImage1 = maskWithinRange(Image1, 1, 350, 1) * Image1
This function clips the range of the data to 1K to 350K. The last "1" passed through the function sets any values outside of this range to be missing. From here, you can pass this masked data object through describe() and you'll see min/max values match the enhancement range returned from the HYDRA chooser.

Code: Select all

print(describe(maskedImage1))
Once you display this maskedImage1 data object, the min/max values will match the values returned from the HYDRA chooser with potentially some minor rounding differences. The enhancement used will be "Gray scale". To match the enhancement used by the HYDRA chooser, you could use setEnhancement() where MOD1 is the layer:

Code: Select all

MOD1.setEnhancement('Inverse Gray Scale')
I'm still looking for a workaround for the missing stripes of data and I'll let you know if I find a solution.

Thanks,
Bob
User avatar
m.yunus.a.m.
Posts: 18
Joined: Wed Aug 18, 2021 3:55 am

Re: How do I customise Default Background Maps?

Post by m.yunus.a.m. »

Hi Bob.

About the issue with dateline. The HYDRA in 1.8 has that problem, but the HYDRA in 1.9beta1 does not have that problem. Perhaps that would help you.

Referring to this post, viewtopic.php?p=7149#p7149. I would like make a suggestion. Perhaps make a custom function that tracks what changes in the UI? For example, Step 1, turn the tracking function on. Step 2, make a change somewhere in the UI. Step 3, turn off the tracking function. Step 4, display what has changes, what functions were called between Step 1 and Step 3. I have no idea if this can be easily implemented or not. So it's just a suggestion. To me, it would be friendlier to user, as it helps user to easily identify what changes in the UI functionality and how to replicate it..

Below is the script I wrote. I would like to have your feedback on it. Technically this is my first Python script, so I think I might not be writing it properly. Basically I have several MOD021km files, and I want to plot the imagery, and the corresponding msl and 10m wind from ERA5 reanalysis (rounded to earlier hour). The data was downloaded earlier based on spatial search on NASA EarthData Search. I took into account of 1 or more imagery in the area, therefore, the script currently checks if there's any granules that are next to each other, and plots both granules, then msl and 10m winds from ERA5.

I have plan to include the possibility of third imagery in the area of interest, but that would be later. That is why there is an empty else. I thought, based on description of 'pass', it would work as place holder, but somehow it doesn't work.

Code: Select all

from datetime import datetime
from ucar.unidata.util import Range

#Grab the current active Display
panel=activeDisplay()

#Open a bundle that has been modified.
openBundle("G:/My Drive/Terra Nova Bay Work in Progress/McIDAS-V/TNB-Ross-Sea.mcv")
#ERA5 file containing mslp, u10, v10, and t2m. The file is hourly.
ERA5Case13='G:/My Drive/Terra Nova Bay Work in Progress/Terra Nova Bay Strong Wind Event/Case13_ERA5_hourly_u10_v10_msl_t2m_20160516-20160529.nc'
#The directory containing the MOD images.
Data_Dir = 'G:/My Drive/Terra Nova Bay Work in Progress/Raw-Data/MODIS-MOD021KM/Terra/20160521'



#Make a list of MODIS imagery in a folder
File_List=sorted(os.listdir(Data_Dir))
#Make a local ADDE entry, not permanent.
Terra20160521 = makeLocalADDEEntry(dataset='TEST1', imageType='MOD021KM', mask=Data_Dir, format='MODS', save=False)  


#Identifier for skipping files. 0 means don't skip. 1 means skips.
skip1=0
#Clear any layers being displayed.
removeAllLayers( )
#for file in range(len(File_List)):
for file in range(len(File_List)):
  
  if skip1==1:
    #If skip was set to 1, do skip, and reset the skip indicator
    print 'Skipping file ',file,' ',File_List[file]
    skip1=0
    continue
  else:
    #pass doesn't work
    print 'Because pass doesn\'t work'
  
  print "Working on file ",file," ", File_List[file]
  #Grab the time details from the MODIS file name.
  file_Time_1 = File_List[file][10:22]
  file_Time_formatted_1 = datetime.strptime(file_Time_1, "%Y%j.%H%M")
  
  
  #Check if it's the last file.
  if file+1==len(File_List):
    file_Time_2=None
    file_Time_formatted_2=None
    continue
  else:
    file_Time_2 = File_List[file+1][10:22]
    file_Time_formatted_2 = datetime.strptime(file_Time_2, "%Y%j.%H%M")
    skip1=1
    #Need to add another in case there are 3 imagery for one area.

  #Try to use dictionary next time.
  Image1=loadADDEImage(server='localhost',localEntry=Terra20160521,\
       time=file_Time_formatted_1.strftime("%X"),day=file_Time_formatted_1.strftime("%Y%j"),\
       band=31,size='ALL',unit='TEMP')
  #To see a summary of the Image1.
  #print(describe(Image1))
  #Mask the image from extremem values.
  maskedImage1 = maskWithinRange(Image1, 1, 350, 1) * Image1
  #Get the new min and max value from masked image.
  #Range1=getMinMax(maskedImage1)
  #To see a summary of masked Image1.
  #print(describe(maskedImage1))
  MOD1=panel.createLayer('Image Display', maskedImage1)
  
  #MOD2.setEnhancement('FROSTWW')
  MOD1.setEnhancement('Inverse Gray shade')
  #Apparently you don't need the line below.
  #MOD1.getGridDisplay().setRangeForColor(Range(Range1))
    
  #If the next file in File_List is the next granule, load it and display it.
  if abs(file_Time_formatted_1.minute-file_Time_formatted_2.minute)==5:
    Image2=loadADDEImage(server='localhost',localEntry=Terra20160521,\
       time=file_Time_formatted_2.strftime("%X"),day=file_Time_formatted_2.strftime("%Y%j"),\
       band=31,size='ALL',unit='TEMP')
    maskedImage2 = maskWithinRange(Image2, 1, 350, 1) * Image2
    MOD2=panel.createLayer('Image Display', maskedImage2)
    #MOD2.setEnhancement('FROSTWW')
    MOD2.setEnhancement('Inverse Gray shade')
    #Need indicator double/triple plotting is needed.
  else:
    #pass doesn't work.
    print 'Because pass doesn\'t work'
  
  #Load ERA5 mslp.
  ERA5MSLP=loadGrid(filename=ERA5Case13,
                  field='msl',latLonBounds=(-55,120,-85,240),time=file_Time_formatted_1.strftime("%Y-%m-%d %H:00:00 UTC"),xStride=1,yStride=1)
  layer1=panel.createLayer('Contour Plan View', ERA5MSLP)
  #Change the display unit into hPa
  layer1.setDisplayUnit('hPa')
  #Setting contour display to be standardised.
  a=layer1.getContourInfo()
  a.interval=4
  a.base=930
  a.max=1030
  layer1.setContourInfo(a)
  #Change the whole shading into red
  layer1.setEnhancement("Red")


  #Loading 10m wind vector.
  ERA5u10=loadGrid(filename=ERA5Case13,
                  field='u10',latLonBounds=(-55,120,-85,240),time=file_Time_formatted_1.strftime("%Y-%m-%d %H:00:00 UTC"),xStride=8,yStride=2)
  ERA5v10=loadGrid(filename=ERA5Case13,
                  field='v10',latLonBounds=(-55,120,-85,240),time=file_Time_formatted_1.strftime("%Y-%m-%d %H:00:00 UTC"),xStride=8,yStride=2)                  


  #Creat display for 10m winds.
  ERA5MSLP10mwind=makeVector(ERA5u10,ERA5v10)
  layer2=panel.createLayer('Vector Plan View', ERA5MSLP10mwind)
  #Change the vecto color to black
  aaa=layer2.getColor()
  layer2.setColor(aaa.BLACK)
  #Standardise the vector range.
  layer2.setFlowRange(Range(-25,25))
  #Save the image.
  panel.captureImage('G:/My Drive/Terra Nova Bay Work in Progress/McIDAS-V/MOD021KM_and_MYD021KM_Imagery/'+file_Time_formatted_1.strftime("%Y%m%d-%H%M")+'.png')
  #Remove all displays for the next iteration.
  removeAllLayers( )
    
print 'Finish'
Yunus.
User avatar
m.yunus.a.m.
Posts: 18
Joined: Wed Aug 18, 2021 3:55 am

Re: How do I customise Default Background Maps?

Post by m.yunus.a.m. »

Hello.

Please help me with putting annotation. Attached is the result of 'title=panel.annotate('2016-05-21 1520 UTC',size=20,color='Black',lat=-67,lon=180,bg='White')'. Is there any way to rotate the annotation? The color scale bar have its orientation correctly according to the display window, but annotate is following the lat lon label orientation.
Attachments
annotate.png
User avatar
m.yunus.a.m.
Posts: 18
Joined: Wed Aug 18, 2021 3:55 am

Re: How do I customise Default Background Maps?

Post by m.yunus.a.m. »

Hello,

Regarding viewtopic.php?p=7363#p7363, and viewtopic.php?p=7377#p7377, it is solved. After contacting IDV support, they suggested to create a new projection and modifying the 'Tangent Lon'. Such as below.
Screen Shot 2021-09-13 at 11.23.15 AM.png
Here's an example of the implementation.
annotate-correct.png
I think Inquiry 2988 and Inquiry 2989 can be closed as the intention behind the requests is covered in this method.

Yunus.
User avatar
bobc
Posts: 987
Joined: Mon Nov 15, 2010 5:57 pm

Re: How do I customise Default Background Maps?

Post by bobc »

Hello,

I'm glad you were able to come up with workarounds for some of your problems. I decided to keep inquiries 2988 and 2989 open as I feel that these are both good suggestions and there may be more general solutions to them other than creating a new projection. I added mention of your workaround to both of these inquiries.

I agree that it would be nice if McIDAS-V could somehow tell the user what the corresponding scripting functions are for User Interface actions. We have received similar requests in the past, and this is written up as inquiry 1197.

As for your script, I would likely need to have all of your data to be able to run the script on my end to give much feedback. From a quick look, it appears that everything in the script is formatted correctly. Were you getting any error messages or were the captured images not what you expected?

Thanks,
Bob
User avatar
m.yunus.a.m.
Posts: 18
Joined: Wed Aug 18, 2021 3:55 am

Re: How do I customise Default Background Maps?

Post by m.yunus.a.m. »

Hi Bob.

The script works fine. I was thinking maybe a more seasoned Python user will find my codes too messy or not up to standard. Or perhaps there's a better option or function to what I have used. I do see some parts that can be simplified by using dictionary declaration.

I do have a different question though. In an attempt to lighten the load on my PC, I tried using File>New Display Window>Misc>2D Map Display. The imagery was plotted above the map lines. It can be seen in the image below. I have also attached the bundle that I created using this "2D Map Display". This doesn't happen with the default 3D display. I'm not sure whether this is working as intended or not, as there is no option to change any vertical position when I use "2D Map Display".
Image being plotted over Map line when using 2D Map Display
Image being plotted over Map line when using 2D Map Display
That being said, I'm going to use the default 3D display for now. There's enough method such as latLonBounds, xStride and yStride to limit the amount of data being loaded.

There is a side issue with using Google Drive (streaming option). The Jython shell sometimes seems to halt. Then the Windows Explorer would stop working. Then any file operation would stop working, even from within other programs. Then I cannot properly shutdown the PC, have to force shutdown the PC. This happens when I'm running the script against the whole directory (which contains at most 17 imagery). I am now working with the files that have been stored locally on my computer instead of cloud (Google Drive). So it's not an issue to me now. But perhaps to other users that uses Google Drive, it will be a problem.

Yunus
Attachments
TNB-Ross-Sea.mcv
(80.48 KiB) Downloaded 85 times
User avatar
bobc
Posts: 987
Joined: Mon Nov 15, 2010 5:57 pm

Re: How do I customise Default Background Maps?

Post by bobc »

Hi Yunus,

I'm replicating the problem of the map lines being hidden behind the data using a 2D map panel if I zoom in close enough to the data. Zoomed out, the map lines appear on top of the data for me. I wrote this up as inquiry 2993. Sometimes this can happen in standard 3D map display panels as well, in which case you could raise the z-positioning of the map lines using the slider at the bottom of the Layer Controls for the map layer. However, this isn't an option with the 2D map display.

Thanks for the information about Google Drive. We've had a couple users in the past ask about using data from Google Drive, but simply URLs pointing to a file on someone's Google Drive (a functionality not currently available in McIDAS-V). I actual never heard of the streaming option, but I just installed it on my end to run some tests.

Thanks again,
Bob
User avatar
m.yunus.a.m.
Posts: 18
Joined: Wed Aug 18, 2021 3:55 am

Re: How do I customise Default Background Maps?

Post by m.yunus.a.m. »

Hello.

I would like to calculate MSLP from WRF output. Since WRF output does not have MSLP by it self, I have to calculate it. I'm referring to http://weather.ou.edu/~criedel/metr2021/plot_wrf.py for the calculation.

Code: Select all

import numpy as np
temps =  nc.variables['T2']
# Convert Surface Pressure to Mean Sea Level Pressure	
stemps = temps[time]+6.5*nc.variables['HGT'][time]/1000.
mslp = nc.variables['PSFC'][time]*np.exp(9.81/(287.0*stemps)*nc.variables['HGT'][time])*0.01 + (6.7 * nc.variables['HGT'][time] / 1000)
In adopting the code to Jython, I referred a work around in https://www.unidata.ucar.edu/support/he ... 08684.html to overcome unit issue.

Code: Select all

def fire(rhscrn,tscrn,u10):
    rhscrn = newUnit(rhscrn,"rhscrn","K")
    tscrn = newUnit(tscrn,"tscrn","K")
    u10 = newUnit(u10,"u10","K")
    answer = 2*exp(1.524-0.0345*rhscrn+0.0338*tscrn+0.0234*u10*3.6)
    return answer
Since Numpy is not available for use, I referred to viewtopic.php?f=26&t=1134, using Numeric to replace Numpy. My attempt is as below.

Code: Select all

import Numeric as np
psfc=loadGrid(filename='',field='PSFC')
temps =  loadGrid(filename='',field='T2')
hgt=loadGrid(filename='',field='HGT')

temps = newUnit(temps,"temps","Pa")
hgt = newUnit(hgt,"hgt","Pa")
psfc = newUnit(psfc,"psfc","Pa")
stemps = temps+6.5*hgt/1000
mslp = psfc*np.exp(9.81/(287.0*stemps)*hgt)*0.01 + (6.7 * hgt / 1000)
It doesn't work however. The issue seems to be with exp in Numeric. It gives the error below

Code: Select all

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: typecode must be in [1silfFdD]
I'm not sure how to proceed with this. I would like help with this.

Yunus
Post Reply