[ncl-talk] help with plotting shapefile polygons
Dave Allured - NOAA Affiliate
dave.allured at noaa.gov
Wed Apr 6 21:39:57 MDT 2022
I found additional problems in one of the NCL script libraries. Here is a
possible solution. Download the attached file
*shapefile_funcs.reversed.temporary.ncl*. Then try this three-line test in
NCL, where *shapefile_utils* is renamed to make it obvious which version is
being used. Specify the *original* shapefile, *not* ogr deprojected. The
order of the load commands is important. Also please ensure that all of
the sidecar files (.prj, .shx, .dbf, etc.) are downloaded into your working
directory, next to the main shapefile:
load "shapefile_funcs.reversed.temporary.ncl"
load "shapefile_utils.reversed.temporary.ncl"
plot_shapefile ("shapefile_name.shp")
I get this PNG file, which looks right:
[image: Perimeter_20180715_SpringCreek_COHUX001313.png]
Vanessa, try this on your system. If it works, then your own program might
also be fixed by adding those two load statements at the top.
On Wed, Apr 6, 2022 at 4:37 PM Rick Brownrigg <brownrig at ucar.edu> wrote:
> I think you need to use the deprojected shapefile that I had sent
> previously.
>
> On Wed, Apr 6, 2022 at 4:35 PM Vanessa Vincente via ncl-talk <
> ncl-talk at mailman.ucar.edu> wrote:
>
>> 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
>>>>>>>>>
>>>>>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/b9b918b2/attachment-0001.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/b9b918b2/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Perimeter_20180715_SpringCreek_COHUX001313.png
Type: image/png
Size: 52298 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/b9b918b2/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: shapefile_funcs.reversed.temporary.ncl
Type: application/octet-stream
Size: 27396 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/b9b918b2/attachment-0001.obj>
More information about the ncl-talk
mailing list