Overlay contour line on a shapefile

[khaled tahoon]

Dear IDL users

I have a scattered datasets of air temperature over Egypt at different meteorological stations.

I would like to plot a contour lines on Egypt' shapefile and exactly inside the Egypt political boundaries.

For this reason, 

firstly, I use a cgScaleVector to create the proper vectors associated with a latitude and longitude coordinate of 51 elements .

secondly,  

the air temperatures are  gridded into a 2D data set of 51 by 51 elements.

   thirdly,

I called Egypt shapefile's inside my own code

my question is

 How can I display the counter lines on Egypt shapefile's and exact only inside the political boundaries.

I wrote the following code and I hope to correct it

------------------------------------------------------------------------------------------------------------------------------------

Pro contour_plot_over_Shapefile

Lat_s=[31.4531,      31.5667,      31.0328,      31.1833,      31.4167,      31.2833,      29.2000,      28.7333,      26.3428,      25.6667,      23.9667,      29.2000,      25.4833,      25.4500,      30.5917, $
           29.8667,     28.2056,      27.9833,      28.9667,      27.1833,         26.9500,22.2, 22.1]
Lon_s=[25.8814,      26.3333,      28.4447,      29.9500,      31.8167,      32.2333,      31.0167,      30.4000,      31.7428,      32.7000,      32.7833,      25.3167,      29.0000,      30.5333,      32.2472, $
            23.9667,      33.6603,      34.3833,      34.6833,      33.8000,     33.9500,31.2,36.1]

Temp=[ 22.1,   25.0,  27.0,   33.0,   18.0,  25.1,   28.0,  15.0, 31.0,  1.0,  24.1,   23.0,  28.0, 35.0,   28.0,   29.1, 38.0,  37.0,   39.0,   38.0,  17.0,34.0,39.0]

; Interpulate the lat and lon data for the shapfile
;------------------------------------------------------------------

lat= Scale_Vector(Findgen(51), 21.99,31.67)     ;          give an array of 51 element
lon = Scale_Vector(Findgen(51), 24.69,36.91)   ;          give an array of 51 element

limit=[min(lat), min(lon), max(lat), max(lon)]
latsubs  = Value_Locate(lat, [limit[0],limit[2]])  ;   Idl print 0     50
lonsubs = Value_Locate(lon, [limit[1],limit[3]]) ;   Idl print 0     50

lat_sub  = lat[latsubs[0]:latsubs[1]]       ; Idl print the lat interpulate over the rang from 0 to 50 and become 51 values
lon_sub  = lon[lonsubs[0]:lonsubs[1]]    ; Idl print the lat interpulate over the rang from 0 to 50 and become 51 values

; --------- gride temp  ( the data )


GRID_INPUT, lon_s, lat_s, temp , xSorted, ySorted, dataSorted

; Initialize the grid parameters.
gridSize = [51,51]
; Use the equation of a straight line and the grid parameters to
; determine the x of the resulting grid.
       slope = (MAX(xSorted) - MIN(xSorted))/(gridSize[0] - 1)
       intercept = MIN(xSorted)
       xGrid = (slope*FINDGEN(gridSize[0])) + intercept

; Use the equation of a straight line and the grid parameters to
; determine the y of the resulting grid.
    slope = (MAX(ySorted) - MIN(ySorted))/(gridSize[1] - 1)
    intercept = MIN(ySorted)
    yGrid = (slope*FINDGEN(gridSize[1])) + intercept

; Grid the data with the Radial Basis Function method.
   grid = GRIDDATA(xSorted, ySorted, dataSorted, $
   DIMENSION = gridSize, METHOD = 'RadialBasisFunction')
;;;;;; print,'kkk', grid
;;;;;;;;;;;; print,  min(grid),  max(grid)
; Levels for contour analysis
nLevels=20
contourLevels = cgConLevels(grid, NLevels=20, MinValue=minValue,  MaxValue=maxValue,Step=step)
tickname=String([contourLevels])

data_sub =grid[lonsubs[0]:lonsubs[1], latsubs[0]:latsubs[1]]
 lat_map=lat(latsubs[0]:latsubs[1],*)  ;--------> y map
 lon_map=lon(lonsubs[0]:lonsubs[1],*)

map_center_lat=(21.99+31.67)/2.0
map_center_lon=( 24.69+36.91)/2.0

cgDisplay, 800, 675
LoadCT, 40, NColors=nlevels, Bottom=1

 cgMap_Set, map_center_lat,map_center_lon,  /CYLINDRICAL , /IsoTropic, Limit=[21.99,24.69,31.67,36.91],    XMARGIN=0, YMARGIN=0,$
Position=[0.05, 0.05, 0.95, 0.95], /NOBORDER                    ; CYLINDRICAL, MERCATOR

shapefile = 'C:\Users\Akram\Desktop\final_shape\shape_file\egy_admbnda_adm0_capmas_itos_20170421.shp'
cgDrawShapes, shapefile, AttrName='ADM0_EN', MapCoord=map, /Projected_XY, Thick=1, $
       FColors=['white'];, Fill = 1;, Colors='white'



  cgContour, grid, xgrid,ygrid,  Levels=contourLevels, Background=cgColor('white'), $
Position=[0.1, 0.1, 0.9, 0.84],    XStyle=1, YStyle=1, C_Colors=IndGen(20)+1, $
ytitle=' Latitude', xtitle='Longitude', /OVERPLOT, Title=' Min. Surfacer Temperature C)'

cgMap_Grid, LatDel = 2.0, LonDel = 2.0, /Box_Axes, Color='charcoal', $
    Map_Structure=map

end

------------------------------------------------------------------------------------

the attach files are Egypt shapefile's

[markus.sc...@gmail.com]

I haven't worked with shapefiles before, so I don't know how to implement it, but it should be possible to set all points in grid outside Egypt to NaN before contouring, this limits the contours to within Egypt. I've used that before in a different context for IDL's contour procedure and function, and it worked. No guarantee for cgContour though. I hope this helps. Has anyone a better idea?

[khaled tahoon]

your idea is a new for me and logic, but you can write how can I set the gridded point outside Egypt to NAN before contouring

[khaled tahoon]

your idea is logic, but how you can set all gridded points outside Egypt shapefile to NAN

if you do something before, could you please write it here

في الاثنين، 17 مايو 2021 في تمام الساعة 7:59:08 م UTC+3، كتب markus.sc...@gmail.com رسالة نصها:

[Ben C]

You should be able to use the lesser known IDLanROI object to determine whether your data being contoured falls within Egypt or not. Specifically the ContainsPoints and ComputeMask methods would be useful, as long as you have all the vertices for the outside boundary of the country. See this Doc page: https://www.l3harrisgeospatial.com/docs/idlanroi.html

Alternatively, you can also do a little trick to have the same effect. Something like this:

• Draw your basemap.
• Overlay all the contoured data for the entire region.
• Overlay an opaque fill of the surrounding countries or bodies of water

The net result is a "window" to view the contours within your ROI of interest. Not sure if it will work with Egypt, but I've done this before for United States. Here's an example. Code below using IDL's built-in function graphics. This may be possible in Coyote graphics, but I am not familiar with them. Hope this helps. 

;IDL SAMPLE CODE FOR CONTOUR WITHIN AN ROI

    ;--------------------------------------------------------------------------------------    

    ;READ VARIABLES FROM FILE

    ;--------------------------------------------------------------------------------------

    no2 = H5_GETDATA(files[i], '/PRODUCT/nitrogendioxide_tropospheric_column')

    time = H5_GETDATA(files[i], '/PRODUCT/time_utc')

    lat = H5_GETDATA(files[i], '/PRODUCT/latitude')

    lon = H5_GETDATA(files[i], '/PRODUCT/longitude')

    no2 = no2 * 1.0e6 

    

    ;Find and use only the data nearby to Colorado and where data is valid (not equal to fill value)

    use=where(no2 lt max(no2) and lat gt 30 and lat lt 50 and lon lt -100 and lon gt -110)

        ;--------------------------------------------------------------------------------------      

        ; DRAW THE BASEMAP

        ;--------------------------------------------------------------------------------------      

        ;SET BOUNDS FOR MAP

          LATMIN=36.5

          LATMAX=41.5

          LONMIN=-109.30

          LONMAX=-101.80

          SCALER = 2.0

        m1 = map('Robinson',LIMIT=[latmin, lonmin, latmax, lonmax],GRID_LONGITUDE=180, GRID_LATITUDE=90,dimensions=[1470*scaler,1000*scaler],$

          title =titletxt,font_style='Bold',MARGIN=0.1)

    

        ;--------------------------------------------------------------------------------------

        ; OVERLAY DATA ONTO MAP VIA CONTOUR

        ;--------------------------------------------------------------------------------------    

        ct = COLORTABLE(72, /reverse)       

        

        c = CONTOUR(no2[use], lon[use],lat[use], /FILL,OVERPLOT=m1, GRID_UNITS='degrees', RGB_TABLE=ct,$

          TITLE='Colorado Tropospheric $NO^{2}$: '+label_time,c_value=indgen(21)*20-10) 

        tit=c.title

        tit.font_size=22

        

        ;--------------------------------------------------------------------------------------    

        ; OVERLAY THE COLORADO COUNTY BOUNDARIES

        ;--------------------------------------------------------------------------------------      

        myshape1 = OBJ_NEW('IDLffShape',co_shapefile)

        ;Get the number of entities so we can parse through them

        myshape1->GetProperty, N_ENTITIES=num_ent

    

        ; Parsing through the entities and only plotting the state of

        ; Colorado

        FOR x=0L, (num_ent-1) DO BEGIN

          ; Get the Attributes for entity x

          attr = myshape1->GetAttributes(x)

          ; See if 'Colorado' is in ATTRIBUTE_1 of the attributes for

          ; entity x

          IF attr.ATTRIBUTE_0 EQ '08' THEN BEGIN

            ; Get entity

            ent = myshape1->GetEntity(x)

            ; Plot entity

            xpoints = (*ent.VERTICES)[0, *]

            ypoints = (*ent.VERTICES)[1, *]

            ; PLOTS, (*ent.VERTICES)[0, *], (*ent.VERTICES)[1, *], COLOR=0

            p = PLOT(xpoints, ypoints,'black',/overplot)

            ; Clean-up of pointers

            myshape1->DestroyEntity, ent

          ENDIF

        ENDFOR

        

        ;--------------------------------------------------------------------------------------    

        ;DRAW STATE LINES NEARBY AND FILL THEM

        ;--------------------------------------------------------------------------------------    

          mc = MAPCONTINENTS(/USA,thick=2,combine=0)

    

            states = ['Wyoming','Utah','Arizona','New Mexico','Oklahoma','Kansas','Nebraska']

            for ss = 0, n_elements(states)-1 do begin

              m1[states[ss]].FILL_COLOR='alice blue'

            endfor