[ncl-talk] help with plotting shapefile polygons
Rick Brownrigg
brownrig at ucar.edu
Wed Apr 6 12:17:26 MDT 2022
Hi Vanessa,
The issue is that the Perimeter_xxxx.shp file's coordinates are in meters
in a Mercator Projection, as seen by the associated Perimeter_xxxx.prj file
(and verified by using NCL to print the x/y coordinates). NCL can only deal
with shapefiles in geographic (lat/lon) coordinates. It's not clear to me
whether the X/Y values are also swapped in this case; I can't tell by
looking at the values.
There are a couple of ways forward: The GDAL package has a command,
ogr2ogr, that in principle, can perform the needed deprojection. Another
possible solution is that NCL has an undocumented interface to the PROJ
cartographic projection library that could do the deproject. Both of these
would take me some experimentation to get it right. I won't be able to give
it a try until later this evening.
Perhaps the best, easiest, and do-it-yourself solution is that if you have
access to GIS tools such as ArcGIS, that could also create a
geographic-coords version of your shapefile.
Rick
On Wed, Apr 6, 2022 at 11:36 AM Vanessa Vincente via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:
> 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.*
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/65c705b9/attachment.html>
More information about the ncl-talk
mailing list