[ncl-talk] Overlay point shapefile on a plot

Karin Meier-Fleischer meier-fleischer at dkrz.de
Mon Jul 3 07:39:56 MDT 2017


Hi Lyndon,

the projection is not the problem I just used the 
cylindrical-equidistant projection to see the global map.

The function gsn_add_shapefile_polymarkers reads the shapefile and scan 
the global file attribute geometry_type. If it doesn't exist it directly 
returns to the main script printing the error message.
There are two ways to go on, you can copy and edit the 
gsn_add_shapefile_polymarkers function or you can draw the polymarkers 
as I did in the modified script using the gsn_add_polymarker function.

Bye, Karin

Am 03.07.17 um 13:24 schrieb Lyndon Mark Olaguera:
> Dear Karin,
>
>
> Thank you for the help. Actually, I was able to plot this successfully 
> by editing the shapefile attribute manually in QGIS. I need to plot it 
> in orthographic projection (Similar to the attached).
>
> Is it possible to specify the "geometry type" within the code (or edit 
> the shapefile attribute?) after reading the shapefile?
>
>
> Many thanks,
>
> Lyndz
>
>
>
> On Mon, Jul 3, 2017 at 8:09 PM, Karin Meier-Fleischer 
> <meier-fleischer at dkrz.de <mailto:meier-fleischer at dkrz.de>> wrote:
>
>     Hi Lyndon,
>
>     the problem is that your shapefile data does not contain point
>     data which is needed by gsn_add_shapefile_polymarkers function
>
>     http://ncl.ucar.edu/Document/Graphics/Interfaces/gsn_add_shapefile_polymarkers.shtml
>     <http://ncl.ucar.edu/Document/Graphics/Interfaces/gsn_add_shapefile_polymarkers.shtml>
>
>     To see what is in your shapefile you can draw the shapefile data
>     using the x- and y-values. I've attached your modified script and
>     the PNG plot.
>
>     Hope this helps,
>     Karin
>
>
>
>     Am 03.07.17 um 04:37 schrieb Lyndon Mark Olaguera:
>>     Dear NCL experts,
>>
>>     I am trying to add point shapefile to a plot using the
>>     gsn_add_shapefile_polymarkers function in NCL.
>>
>>     To view the contents of the shapefile I used the shape_utils.ncl
>>     from the following link:
>>
>>     *_https://www.ncl.ucar.edu/Applications/shapefiles.shtml
>>     <https://www.ncl.ucar.edu/Applications/shapefiles.shtml>
>>     _*
>>
>>     I tried to plot the shapefile in GIS and it works fine but I
>>     still got an error that the geometry attribute should be a point
>>     when I plot it in NCL. I am not sure why in the filedump below,
>>     the geometry type is "unknown".
>>
>>
>>     I attached my script and the files that I am plotting.
>>
>>
>>     Variable: a
>>     Type: file
>>     File path:may_hgt-hgt_1996-2012_850hPa.nc
>>     Number of global attributes:0
>>     Number of dimensions:2
>>     Number of variables:6
>>     (0)======================================================================
>>     (0)Filename: "sig_regpval_1996-2012.shp"
>>     (0)  Geometry type: unknown
>>     (0)  # of features: 2906
>>     (0)  Min/max lat:    -80.00/  82.50
>>     (0)  Min/max lon:      0.00/ 360.00
>>     (0)  Variable names and their types:
>>     (0)      geometry : integer
>>     (0)      segments : integer
>>     (0)      x : double
>>     (0)      y : double
>>     (0)      z : double
>>     (0)      CREATED_BY : string
>>     (0)      AUTHOR : string
>>     (0)      TYPE : string
>>     (0)      DESC : string
>>     (0)      LONGITUDE : double
>>     (0)      LATITUDE : double
>>     (0)      GRID_VALUE : double
>>     (0)======================================================================
>>     (0)Error: gsn_add_shapefile_polymarkers: geometry_type attribute
>>     must be 'point'
>>     (0)      No shapefile information will be added.
>>
>>
>>     I'll appreciate any suggestion on how I can plot this correctly.
>>
>>     Thank you in advance.
>>
>>     Lyndz
>>
>>
>>
>>     _______________________________________________
>>     ncl-talk mailing list
>>     ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
>>     List instructions, subscriber options, unsubscribe:
>>     http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>     <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
>
>     -- 
>     Dipl. Geophys. Karin Meier-Fleischer
>     Visualization, NCL, CDO
>     Application Support
>
>     Deutsches Klimarechenzentrum GmbH (DKRZ)
>     Bundesstrasse 45a - D20146 Hamburg - Germany
>
>     Phone:+49 (0)40 460094 126 <tel:+49%2040%20460094126>
>     Fax:+49 (0)40 460094 270 <tel:+49%2040%20460094270>
>     E-Mail:meier-fleischer at dkrz.de <mailto:meier-fleischer at dkrz.de>
>     URL:www.dkrz.de <http://www.dkrz.de>
>
>     Geschäftsführer: Prof. Dr. Thomas Ludwig
>     Sitz der Gesellschaft: Hamburg
>     Amtsgericht Hamburg HRB 39784
>
>     _______________________________________________ ncl-talk mailing
>     list ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu> List
>     instructions, subscriber options, unsubscribe:
>     http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>     <http://mailman.ucar.edu/mailman/listinfo/ncl-talk> 
>
-- 
Dipl. Geophys. Karin Meier-Fleischer
Visualization, NCL, CDO
Application Support

Deutsches Klimarechenzentrum GmbH (DKRZ)
Bundesstrasse 45a - D20146 Hamburg - Germany

Phone:    +49 (0)40 460094 126
Fax:      +49 (0)40 460094 270
E-Mail:   meier-fleischer at dkrz.de
URL:      www.dkrz.de

Geschäftsführer: Prof. Dr. Thomas Ludwig
Sitz der Gesellschaft: Hamburg
Amtsgericht Hamburg HRB 39784
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170703/0a503270/attachment-0001.html 
-------------- next part --------------
begin

  a = addfile("may_hgt-hgt_1996-2012_850hPa.nc","r")
  
;************************************************
; read in data
;************************************************
  printVarSummary(a)
  t = a->regslope(:,:)   
                               
;************************************************
; create plot
;************************************************
  wks = gsn_open_wks("png","proj")                ; send graphics to PNG file
       
  res                   = True      
  res at gsnFrame          = False                   ;-- don't advance the frame
        
  res at cnLevelSelectionMode = "ManualLevels"      
  res at cnMinLevelValF    = -5.00      
  res at cnMaxLevelValF    = 5.00      
  res at cnLevelSpacingF   = 0.5      
      
  res at mpProjection      = "Orthographic"          ; choose projection
  res at mpCenterLonF      = 122.                    ; choose center lon
  res at mpCenterLatF      = 45.                     ; choose center lat
  res at mpPerimOn         = False                   ; turn off box around plot
  res at mpFillOn          = False      
  res at mpDataBaseVersion = "MediumRes"      
  res at mpOutlineOn       = True      
      
  res at cnFillOn          = True                    ; color plot desired
  res at cnFillPalette     = "BlueDarkRed18"           ; select color map
  res at cnLineLabelsOn    = False                   ; turn off contour line labels
  res at cnLinesOn         = False                   ; turn off contour lines
  
  res at mpGridAndLimbOn        = True               ; turn on lat/lon grid lines
  res at mpGridMaskMode         = "MaskNotOcean"     ; don't draw over land or
                                                  ; inland water bodies

;  res at lbLabelFontHeightF  = 0.015                 ; label bar font height
 
  res at tiMainString       = "Test Orthographic Projection"   ; add a title
;  res at tiMainFontHeightF  = .018                   ; font height

  plot = gsn_csm_contour_map(wks,t,res)

;;************************************************************
;;Add shapefiles
;;************************************************************

  shpname  = "sig_regpval_1996-2012.shp"
  sf       =  addfile(shpname,"r")
  
  vNames   =  getfilevarnames(sf)                  ;-- variable names
  x        =  sf->x
  y        =  sf->y
  npoints  =  dimsizes(x)                          ;-- number of markers

;-- polymarker resources
  pres                  =  True
  pres at gsMarkerIndex    =  16         
  pres at gsMarkerColor    = "black"
  pres at gsMarkerOpacityF =  0.3                     ;-- make the marker more transparent
  pres at gsMarkerSizeF    =  0.005

;-- draw all polymarkers
  do i=0,npoints-1
     id = unique_string("plid")
     plot@$id$ = gsn_add_polymarker(wks,plot,x(i),y(i),pres)
  end do
  
  draw(plot)
  frame(wks)
end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: proj.png
Type: image/png
Size: 979127 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170703/0a503270/attachment-0001.png 


More information about the ncl-talk mailing list