[ncl-talk] help with plotting shapefile polygons

[Vanessa Vincente]

Hello,

I am having trouble plotting polygon data from a shapefile of a burn scar
perimeter over a map of Colorado with counties.  I am using the example
shapefiles_15.ncl from the ncl page
<https://www.ncl.ucar.edu/Applications/shapefiles.shtml>. I am able to
successfully plot the state of Colorado and its counties, but not the burn
scar perimeter.  Burn scar shapefile can be found here
<https://ftp.wildfire.gov/public/incident_specific_data/rocky_mtn/2018/SpringCreek/IR/20180715/PerimeterFileProvidedByIncident/>.
 Can you please help? Thanks in advance!

Below is a copy of my ncl script:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "./shapefile_utils.ncl"

;----------------------------------------------------------------------
; This function creates a cylindrical equidistant map of Colorado.
;----------------------------------------------------------------------

function create_colorado_map(wks,res,draw_ncl_outlines)
local a, mpres
begin
  mpres                       = res

  mpres at gsnMaximize           = True
  mpres at gsnPaperOrientation   = "portrait"

  mpres at gsnDraw               = False
  mpres at gsnFrame              = False

  mpres at mpFillOn              = False

;---Turn on fancier tickmark labels.
  mpres at pmTickMarkDisplayMode = "Always"
  mpres at tmXBLabelFontHeightF  = 0.008      ; smaller tickmark labels

;---Zoom in on area of interest
  mpres at mpLimitMode           = "LatLon"
  mpres at mpMinLatF             =   37
  mpres at mpMaxLatF             =   41
  mpres at mpMinLonF             = -109.05
  mpres at mpMaxLonF             = -102.05
  mpres at mpFillOn              = False

  if(draw_ncl_outlines) then
    mpres at mpOutlineOn           = True
    mpres at mpOutlineBoundarySets    = "AllBoundaries"
    mpres at mpDataBaseVersion        = "MediumRes"
    mpres at mpDataSetName            = "Earth..4"      ; U.S. counties
  else
    mpres at mpOutlineOn              = False
  end if

;---Create map.
  map = gsn_csm_map(wks,mpres)

  return(map)
end

;--------------------------------------------------
; Main code
;--------------------------------------------------
begin
  wtype          = "png"              ; send graphics to PNG file
  wtype at wkWidth  = 2000
  wtype at wkHeight = 2000
  wks = gsn_open_wks(wtype,"shapefiles")

  ncl_version = get_ncl_version()

;---Create two maps of Colorado
  res                   = True
  res at tiMainFontHeightF = 0.015

  res at tiMainString  = "Colorado counties - shapefile"
  map_shp = create_colorado_map(wks,res,False)

  res at tiMainString = "Colorado counties - NCL (version " + ncl_version + ")"
  map_ncl  = create_colorado_map(wks,res,True)

;---Add shapefiles to one of the maps
  lnres             = True
  lnres at gsLineColor = "red"

  plot_line_shp =
gsn_add_shapefile_polylines(wks,map_shp,"./usa/gadm40_USA_2.shp",lnres)

;---Add a slightly transparent marker to both maps to show location of
county updates
  mkres                  = True
  mkres at gsMarkerIndex    = 16        ; filled dot
  mkres at gsMarkerOpacityF = 0.5       ; make the marker half transparent
  mkres at gsMarkerColor    = "red"
  mkres at gsMarkerSizeF    = 50.

  counties_lat_center =   39.8
  counties_lon_center = -104.9

  plot_marker_ncl =
gsn_add_polymarker(wks,map_ncl,counties_lon_center,counties_lat_center,mkres)
  plot_marker_shp =
gsn_add_polymarker(wks,map_shp,counties_lon_center,counties_lat_center,mkres)

print("Adding polygons to ncl map and shp map...")

filename = "Perimeter_20180715_SpringCreek_COHUX001313.shp"

 ;---Set some options for the polygons
   bscar = True
   bscar at gsLineColor = "Blue"
   bscar at gsLineThicknessF = 3.0  ; 3x thickness

id1 = gsn_add_shapefile_polygons(wks,map_ncl,filename,bscar)
id2 = gsn_add_shapefile_polygons(wks,map_shp,filename,bscar)

;---Panel both plots. Markers, lines, and polygons will be drawn.
  pres             = True
  pres at gsnMaximize = True
  gsn_panel(wks,(/map_ncl,map_shp/),(/1,2/),pres)

end



-- 
Vanessa Vincente
Associate Scientist
The COMET® Program
University Corporation for Atmospheric Research

*My working hours may differ from your own.  Please do not feel obligated
to reply outside of standard business hours.*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220405/648eb316/attachment.html>