[ncl-talk] help with plotting shapefile polygons
Vanessa Vincente
vincente at ucar.edu
Wed Apr 6 15:50:25 MDT 2022
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.*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220406/eb6d2f5f/attachment.html>
More information about the ncl-talk
mailing list