[ncl-talk] help with plotting shapefile polygons
Rick Brownrigg
brownrig at ucar.edu
Wed Apr 6 16:37:14 MDT 2022
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
>>>>>>>>
>>>>>>>>
>>>
>>> --
>>> 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.*
> _______________________________________________
> 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/e2d0499c/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/e2d0499c/attachment.png>
More information about the ncl-talk
mailing list