[ncl-talk] help with plotting shapefile polygons
Vanessa Vincente
vincente at ucar.edu
Wed Apr 6 11:35:30 MDT 2022
Thank you for sharing this workaround. I have replaced the
shapefile_utils.ncl with the one provided, and ran the code, but it is
still not displaying the burn scar perimeter polygon. I am running the ncl
code in the directory that houses all the shapefiles and the new
shapefile_utils.ncl.
I am using version 6.6.2 and the source of the NCL installation was conda
4.11.0.
Below is a copy of my code again for reference.
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
On Tue, Apr 5, 2022 at 4:39 PM Dave Allured - NOAA Affiliate <
dave.allured at noaa.gov> wrote:
> You may have encountered a known bug whereby the X and Y coordinates from
> a shapefile are swapped with each other. We think this was caused by a
> change in an underlying library. Please see this report for details and a
> possible easy workaround. Let us know if that works.
>
> https://github.com/NCAR/ncl/issues/176
>
> What version of NCL are you using, and what was the source of your NCL
> installation?
>
>
> On Tue, Apr 5, 2022 at 4:29 PM Vanessa Vincente via ncl-talk <
> ncl-talk at mailman.ucar.edu> wrote:
>
>> 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.*
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at mailman.ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
--
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/20220406/56c7336e/attachment.html>
More information about the ncl-talk
mailing list