Does CONTOUR /OVERPLOT /FILL work on MAP_SET?

Help
2014-01-16
2014-01-24
  • Hyo-Kyung Lee
    Hyo-Kyung Lee
    2014-01-16

    Hi,

    I have a question about CONTOUR /FILL /OVERPLOT.

    The following code works beautifully and generates color contour as I expected:

    url='http://eosdap.hdfgroup.uiuc.edu:8080/opendap/data/NASAFILES/hdf4/AIRS.2002.08.01.L3.RetStd_H031.v4.0.21
    .0.G06104133732.hdf'
    fid = NCDF_OPEN(url)
    vid = NCDF_VARID(fid, 'Temperature_MW_A')
    NCDF_ATTGET, fid, vid, 'long_name', long_name
    PRINT, long_name
    NCDF_ATTGET, fid, vid, '_FillValue', fillvalue
    PRINT, fillvalue
    NCDF_VARGET, fid, vid, data
    
    ; Subset data.
    data2D = data[*,*,11]
    data2D=Reform(data2D)
    
    ; Process fill values, convert data that are equal to fillvalue to NaN
    dataf=float(data2D)
    fillvaluef = float(fillvalue(0))
    idx=WHERE(dataf eq fillvaluef, cnt)
    IF cnt GT 0 THEN dataf[idx] = !Values.F_NAN
    
    vid = NCDF_VARID(fid, 'Longitude')
    NCDF_VARGET, fid, vid, lon
    
    vid = NCDF_VARID(fid, 'Latitude')
    NCDF_VARGET, fid, vid, lat
    
    levels=254
    DEVICE
    LOADCT, 33, NCOLORS=levels, BOTTOM=1
    CONTOUR, dataf, lon, lat, /FILL, NLEVELS=levels, POSITION=[0.05, 0.06, 0.82, 0.80]
    
    units = 'K'
    XYOUTS, 0.5, 0.9, /NORMAL, url , CHARSIZE=1.3, ALIGNMENT=0.5
    XYOUTS, 0.05, 0.86, /NORMAL, 'FIELD:Temperature_MW_A at TempPrsLvls=11', $
        CHARSIZE=1.25,  ALIGNMENT=0.0
    XYOUTS, 0.94, 0.86, /NORMAL, 'UNIT:' + units, $ 
        CHARSIZE=1.25,  ALIGNMENT=1.0
    

    However, if I add MAP_SET and /OVERPLOT on CONTOUR for the same code, it gives an error.

    PLPLOT ERROR, ABORTING OPERATION
    plfshade: shade_max must exceed shade_min, aborting operation

    At the end, GDL shows only box grid and map continent.
    Doesn't CONTOUR /FILL work with MAP_SET yet in GDL 0.9.4 with PLplot 5.9.11?

    url='http://eosdap.hdfgroup.uiuc.edu:8080/opendap/data/NASAFILES/hdf4/AIRS.2002.08.01.L3.RetStd_H031.v4.0.21
    .0.G06104133732.hdf'
    fid = NCDF_OPEN(url)
    vid = NCDF_VARID(fid, 'Temperature_MW_A')
    NCDF_ATTGET, fid, vid, 'long_name', long_name
    PRINT, long_name
    NCDF_ATTGET, fid, vid, '_FillValue', fillvalue
    PRINT, fillvalue
    NCDF_VARGET, fid, vid, data
    
    ; Subset data.
    data2D = data[*,*,11]
    data2D=Reform(data2D)
    
    ; Process fill values, convert data that are equal to fillvalue to NaN
    dataf=float(data2D)
    fillvaluef = float(fillvalue(0))
    idx=WHERE(dataf eq fillvaluef, cnt)
    IF cnt GT 0 THEN dataf[idx] = !Values.F_NAN
    
    vid = NCDF_VARID(fid, 'Longitude')
    NCDF_VARGET, fid, vid, lon
    
    vid = NCDF_VARID(fid, 'Latitude')
    NCDF_VARGET, fid, vid, lat
    
    levels=254
    DEVICE
    LOADCT, 33, NCOLORS=levels, BOTTOM=1
    MAP_SET, /GRID, /CONTINENT, POSITION=[0.05, 0.06, 0.82, 0.80]
    CONTOUR, dataf, lon, lat, /OVERPLOT, /FILL, NLEVELS=levels, POSITION=[0.05, 0.06, 0.82, 0.80]
    units = 'K'
    XYOUTS, 0.5, 0.9, /NORMAL, url , CHARSIZE=1.3, ALIGNMENT=0.5
    XYOUTS, 0.05, 0.86, /NORMAL, 'FIELD:Temperature_MW_A at TempPrsLvls=11', $
        CHARSIZE=1.25,  ALIGNMENT=0.0
    XYOUTS, 0.94, 0.86, /NORMAL, 'UNIT:' + units, $ 
        CHARSIZE=1.25,  ALIGNMENT=1.0
    MAP_GRID, /BOX_AXES, COLOR=WHITE, /HIRES
    MAP_CONTINENTS, COLOR=WHITE, /HIRES
    
     
  • Alain C.
    Alain C.
    2014-01-16

    Hi

    thanks for the feedback, we learned a lot !
    Do you really succeed to run in GDL:

    url='http://eosdap.hdfgroup.uiuc.edu:8080/opendap/data/NASAFILES/hdf4/AIRS.2002.08.01.L3.RetStd_H031.v4.0.21.0.G06104133732.hdf'
    fid = NCDF_OPEN(url)

    If yes, please detailed the OS + netcdf version !
    I tried on various recent Linux machines, without success up to now.
    I read it is related to opendap, now available in netcdf 4.1+
    Never tried before.

    (I am trying to reproduce running you code)

    Alain

    PS: I am working on NetCDF-4 functionalities in GDL
    since a week, I would really appreciate feedbacks
    to test and improve the code (e.g. I have no example
    for NCDF_UNLIMDIMSINQ)

     
  • Alain C.
    Alain C.
    2014-01-17

    OK, I understood: I broke few days ago this DAP in netCDF trying to improve
    the way NCDF_OPEN() checks files !

    --> I will correct it ASAP (today)

    --> I will add few extra tests cases !

    thanks

    Alain

    PS: after that I will give a look to your question !

     
  • giloo
    giloo
    2014-01-17

    Hi,
    Thanks for this nice example that will be useful for me to test CONTOUR in the future.

    I can reproduce your bug.
    Actually several plotting routines in 0.9.4 behave better than in 0.9.3 but... are broken when using projections. We'll have to investigate/repair, this is known. Apparently CONTOUR not only ignores projection but destroys the projection information.

    You could try to CONTOUR but with map already projected.
    I have successfully done (but map_patch is terribly long with GDL):
    MAP_SET, POSITION=[0.05, 0.06, 0.82, 0.80]
    patched=map_patch(dataf,lon,lat)
    CONTOUR, patched, POSITION=[0.05, 0.06, 0.82, 0.80],nlev=24 ; kills PROJ ^(
    MAP_SET, POSITION=[0.05, 0.06, 0.82, 0.80],/noerase ; restore proj
    MAP_CONTINENTS, COLOR=WHITE, /HIRES

    I noticed also that MAP_PROJ_FORWARD hangs GDL. Another projection-related bug.

    For Alain C.: I too succeed in loading the data with ncdf_open (version { 0.9.4 CVS Jan 5 2014 1388876400} on an up to date mageia3 distro).

     
  • giloo
    giloo
    2014-01-23

    Hi,
    Actually I can reproduce your example. it was not using map_patch, i was wrong.
    The problem is, since Contour 'ignores' (hopefully temporarily) the existence of a projection, you must map exactly the coordinates of the 2D contour image (vectors lon,lat) with the contours of continents.
    just do:

    MAP_SET, POSITION=[0.05, 0.06, 0.82, 0.8] ; set projection
    CONTOUR, dataf, lon, lat, /FILL, NLEVELS=levels, POSITION=[0.05, 0.06, 0.82, 0.80],/xst,/yst,xtickint=45,ytickint=30 ; contour kills projection info
    MAP_SET, POSITION=[0.05, 0.06, 0.82, 0.8],/noerase ; restore projection info
    MAP_CONTINENTS, COLOR=WHITE ; draw continents
    

    I used loadct,39 for the attached plot, since loadct,33 plots the text in black with my configuration (?).

     
    Attachments
    • Hyo-Kyung Lee
      Hyo-Kyung Lee
      2014-01-24

      Thanks, giloo! This is awesome!

      I could get the same nice plot over map so I updated the GDL webpage with a new screenshot and code:

      http://hdfeos.org/software/gdl.php

      I put an acknowledgement in the complete sample code:

      http://hdfeos.org/examples/code/gdl/opendap.gdl

      If you'd like to see your real name and organization in the acknowledgement section or other acknowledgement statement, please feel free to contact me at eoshelp@hdfgroup.org. I really appreciate your help and support.

      p.s. I hope GDL can package map support as a built-in feature. Building from the source was too difficult and error-prone.