[ncl-talk] help with plotting shapefile polygons
Rick Brownrigg
brownrig at ucar.edu
Wed Apr 6 12:38:48 MDT 2022
Actually, it turned out to be quite easy to deproject the file with ogr2ogr:
ogr2ogr -t_srs EPSG:4326
Perimeter_20180715_SpringCreek_COHUX001313_latlon.shp
Perimeter_20180715_SpringCreek_COHUX001313.shp
The resultant coords are definitely in lat/lon and look like values for
Colorado. I *assume* ogr2gr used the information in the *.prj file, to
properly set up the requisite deprojection parameters -- check that the
outline is spatially registered correctly. Also, the swapped X/Y *may*
still be an issue, I can't tell.
Rick
On Wed, Apr 6, 2022 at 12:17 PM Rick Brownrigg <brownrig at ucar.edu> wrote:
> 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/e32ea256/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PerimLL.zip
Type: application/zip
Size: 163059 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/e32ea256/attachment.zip>
More information about the ncl-talk
mailing list