[ncl-talk] help with plotting shapefile polygons

Vanessa Vincente vincente at ucar.edu
Wed Apr 6 16:35:17 MDT 2022


Thanks Dave.  I ran the diagnostic with the replacement
shapefile_utils.ncl, and checked the PNG it produced. Unfortunately,
nothing was plotted in the PNG, see below. Does this mean that NCL is
having trouble understanding the shapefile?
[image: image.png]

On Wed, Apr 6, 2022 at 4:13 PM Dave Allured - NOAA Affiliate <
dave.allured at noaa.gov> wrote:

> Those are geographic coordinates, approximately correct for the Spring
> Fire; and they are swapped.  Therefore, use the *replacement
> shapefile_utils.ncl*; and *do not* convert any farther with ogr2ogr.
>
> Use this diagnostic to make a simple check plot for the whole perimeter.
> This is from the shapefiles example page, under "Shapefile plotting and
> masking routines".  This will find out whether NCL can fully understand
> this shapefile and place it on a map.  This is silent, but it should
> immediately write a PNG file in your working directory.  Open the PNG file
> with any browser or graphics viewer.
>
>     load "./shapefile_utils.ncl"           (the replacement version)
>     plot_shapefile("shapefile_name.shp")
>
>
> On Wed, Apr 6, 2022 at 3:50 PM Vanessa Vincente <vincente at ucar.edu> wrote:
>
>> Dave and Rick,
>>
>> Thank you both for your help.  I have conducted the test locally, and get
>> the following:
>>
>> Variable: x
>>
>> Type: double
>>
>> Total Size: 110992 bytes
>>
>>             13874 values
>>
>> Number of Dimensions: 1
>>
>> Dimensions and sizes: [num_points | 13874]
>>
>> Coordinates:
>>
>>             num_points: not a coordinate variable
>>
>> Number of Attributes: 0
>>
>>
>> (0) Print MinMax for x...
>>
>> (0) min=37.39254516610051   max=37.69200009599813
>>
>> (0) Print MinMax for y...
>>
>> (0) min=-105.2910941200151   max=-105.0002034698655
>>
>>
>> I did this test using both the original shapefile_utils.ncl and
>> replacement shapefile_untils.ncl, and got the same output. It looks like
>> latitudes are x and latitudes are y. Does this mean I can use the original
>> shapefile_utils.ncl, but I still need to deproject the Perimeter_...shp
>> file using ogr2ogr?  If so, can you please provide a bit more guidance on
>> how I can deproject?
>>
>> Thanks,
>> Vanessa
>>
>>
>>
>> On Wed, Apr 6, 2022 at 1:02 PM Dave Allured - NOAA Affiliate <
>> dave.allured at noaa.gov> wrote:
>>
>>> Rick, thank you for discovering the projection problem.  Here is a
>>> simple way to determine whether coordinates are swapped.  Remember, the
>>> swapping bug is a characteristic of the local NCL installation, not the
>>> shapefile.  So Vanessa will need to do this test locally.  Exclude all
>>> outside influences, and verify the reprojected shapefile by NCL direct
>>> access:
>>>
>>>     f = addfile (shapefile_name, "r")
>>>     printFileVarSummary(f,"x")
>>>     printMinMax (f->x, 0)
>>>     printMinMax (f->y, 0)
>>>
>>> If geographic projection, then longitudes are x, and latitudes are y.
>>> This will immediately tell you whether NCL is reading normal or swapped
>>> coordinates, also whether the shapefile is even valid or readable.  It will
>>> also determine metric vs. geographic coordinates, if you understand what
>>> ranges to expect.
>>>
>>> If the coordinates are normal from this test, then use the original
>>> shapefile_utils.ncl.  If the coordinates are swapped, then use the
>>> replacement shapefile_utils.ncl.
>>>
>>>
>>> On Wed, Apr 6, 2022 at 12:39 PM Rick Brownrigg <brownrig at ucar.edu>
>>> wrote:
>>>
>>>> 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
>>>>>>>> _______________________________________________
>>>>>>>> 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.*
>>
>

-- 
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/5adbac2e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 60808 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/5adbac2e/attachment.png>


More information about the ncl-talk mailing list