[ncl-talk] Plotting Lat/Lon Points on a map
Mary Haley
haley at ucar.edu
Thu Jun 1 09:37:57 MDT 2017
Hi Dave,
We have a ticket (NCL-2207) for modernizing the map outlines in NCL, and we
hope to do this through the use of public domain sources, like GADM or
similar. I just made this ticket a higher priority, and it recently was put
in our road map as a high priority.
I agree that the use of shapefiles is the best route for now. We have many
examples at:
http://www.ncl.ucar.edu/Applications/shapefiles.shtml
and a template for adding shapefile outlines to an existing plot:
http://www.ncl.ucar.edu/Applications/Templates/#AddingShapefileOutlinesTemplates
--Mary
On Tue, May 30, 2017 at 7:32 PM, Dave Allured - NOAA Affiliate <
dave.allured at noaa.gov> wrote:
> Tim and NCL support,
>
> State and county lines in the current NCL database are not accurate for
> high resolution maps. Please see the attached figure. Near Tim's
> reference point, NCL's built-in NC/VA state line is about 1 kilometer north
> of the same line in the GADM United States shape file. This offset is
> large enough to move the appearance of the reference point into the wrong
> state. It appears that the GADM shape file is much more accurate than the
> NCL database.
>
> Can the NCL built-in database be updated with more accurate state and
> county lines?
>
> Tim, for an interim fix I recommend disabling the built-in state
> outlines. Instead, plot state and county outlines using an accurate shape
> file such as the GADM U.S. shape file. See the attached NCL script for an
> example. I recommend against using any manual offsets, because those can
> lead to abundant hidden mistakes.
>
> In addition, I also found discrepancies in state and county lines, between
> the GADM U.S. shape file and Google Earth. The state line discrepancy at
> this Virginia location is about 100 meters, or in general about an "order
> of magnitude" less than the same discrepancy with NCL.
>
> --Dave
>
>
> On Sun, May 28, 2017 at 10:13 AM, Tim Melino <melino33 at gmail.com> wrote:
>
>> Thanks Barry and Dennis for your help on this!
>>
>> Dennis, I added your my other question was why is the point in NCL
>> plotted south of the North Carolina border with NCL? The actual location is
>> in Virginia, I noticed this was occurring in the image that you sent me
>> too.
>>
>> [image: Inline image 1]
>>
>> On Sun, May 28, 2017 at 12:00 PM, Barry Lynn <barry.h.lynn at gmail.com>
>> wrote:
>>
>>> Hi Dennis:
>>>
>>> Thank you for helping. Sorry I was unavailable to (try to) do more.
>>>
>>> Barry
>>>
>>> On Sun, May 28, 2017 at 6:25 PM, Dennis Shea <shea at ucar.edu> wrote:
>>>
>>>> See attached. Please look carefully at the output produced by
>>>> printVarSummary and printMinMax
>>>>
>>>> Note:
>>>>
>>>> [1] GRIB is a moving target. 6.3.0 did not have the look-up tabl
>>>> builtin while NCL 6.4.0 does.
>>>>
>>>> [2] the grib file has 'no report' coded as as -1 ... the script sets
>>>> those values to _FillValue
>>>>
>>>> [3] The grid is *rectilinear* :
>>>> In netCDF speak, a one dimensional variable where the variable name is
>>>> the same as the variable dimension name is a coordinate variable.
>>>> lat_0(lat_0) and lon_0(lon_0)
>>>>
>>>> Good Luck
>>>>
>>>>
>>>>
>>>> On Sun, May 28, 2017 at 10:21 AM, Barry Lynn <barry.h.lynn at gmail.com>
>>>> wrote:
>>>>
>>>>> Tim:
>>>>>
>>>>> If you can wait a day or two, I suggest you write to:
>>>>>
>>>>> Carrie.Langston at noaa.gov
>>>>>
>>>>> and ask Him/Her what is the format of the gridded data and if they
>>>>> have an NCL program to read it.
>>>>>
>>>>> It's looks like very useful data once you "crack the code."
>>>>>
>>>>> Barry
>>>>>
>>>>> On Sun, May 28, 2017 at 5:09 PM, Tim Melino <melino33 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Barry,
>>>>>>
>>>>>> The source of the data is here:
>>>>>>
>>>>>> http://www.nssl.noaa.gov/projects/mrms/operational/tables.php
>>>>>>
>>>>>>
>>>>>>
>>>>>> Tim
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Sun, May 28, 2017 at 10:03 AM, Barry Lynn <barry.h.lynn at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> Tim:
>>>>>>>
>>>>>>> If you write out the latitude and longitude information is is
>>>>>>> rectilinear?
>>>>>>>
>>>>>>> When trying to sort this out, you should first check if the lat/lon
>>>>>>> can be plotted by a single coordinate or two. If two, then there is a
>>>>>>> chance that your grid is not rectilinear, and you will have to account for
>>>>>>> the distortion in the grid.
>>>>>>>
>>>>>>> If there is distortion, then you need to account for it (e.g.,
>>>>>>> Lambert, curvilinear.
>>>>>>>
>>>>>>> I can only offer advice at the moment...
>>>>>>>
>>>>>>> What is the data source for the data?
>>>>>>>
>>>>>>> Barry
>>>>>>>
>>>>>>> On Sun, May 28, 2017 at 4:52 PM, Tim Melino <melino33 at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hey Barry,
>>>>>>>>
>>>>>>>> The only information about the projection was: *grid_type :
>>>>>>>> Latitude/Longitude. *I attached the file!
>>>>>>>>
>>>>>>>> Tim
>>>>>>>>
>>>>>>>> On Sun, May 28, 2017 at 9:00 AM, Barry Lynn <barry.h.lynn at gmail.com
>>>>>>>> > wrote:
>>>>>>>>
>>>>>>>>> Tim:
>>>>>>>>>
>>>>>>>>> Please remind me what the name of the file is.
>>>>>>>>>
>>>>>>>>> Also, was there any information at the top of the printing about
>>>>>>>>> the projection?
>>>>>>>>>
>>>>>>>>> Yes, the markers would (could) look off.
>>>>>>>>>
>>>>>>>>> On Sun, May 28, 2017 at 3:32 PM, Tim Melino <melino33 at gmail.com>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Hi Barry,
>>>>>>>>>>
>>>>>>>>>> Here is the information when printing (a). Would the map
>>>>>>>>>> projection also cause a marker to be printed in the wrong location?
>>>>>>>>>>
>>>>>>>>>> Variable: a
>>>>>>>>>> Type: file
>>>>>>>>>> filename: test
>>>>>>>>>> path: ./test.grib2
>>>>>>>>>> file global attributes:
>>>>>>>>>> dimensions:
>>>>>>>>>> lat_0 = 30
>>>>>>>>>> lon_0 = 40
>>>>>>>>>> variables:
>>>>>>>>>> float VAR_209_3_34_P0_L102_GLL0 ( lat_0, lon_0 )
>>>>>>>>>> center : US NOAA Office of Oceanic and Atmospheric
>>>>>>>>>> Research
>>>>>>>>>> production_status : Research products
>>>>>>>>>> long_name : unknown variable name
>>>>>>>>>> units : unknown
>>>>>>>>>> _FillValue : 1e+20
>>>>>>>>>> grid_type : Latitude/longitude
>>>>>>>>>> parameter_discipline_and_category : 209, 3
>>>>>>>>>> parameter_template_discipline_category_number :
>>>>>>>>>> ( 0, 209, 3, 34 )
>>>>>>>>>> level_type : Specific altitude above mean sea level (m)
>>>>>>>>>> level : 500
>>>>>>>>>> forecast_time : 0
>>>>>>>>>> forecast_time_units : minutes
>>>>>>>>>> initial_time : 07/09/2016 (04:46)
>>>>>>>>>>
>>>>>>>>>> float lat_0 ( lat_0 )
>>>>>>>>>> long_name : latitude
>>>>>>>>>> grid_type : Latitude/Longitude
>>>>>>>>>> units : degrees_north
>>>>>>>>>> Dj : 0.01
>>>>>>>>>> Di : 0.009999676
>>>>>>>>>> Lo2 : 280.735
>>>>>>>>>> La2 : 36.695
>>>>>>>>>> Lo1 : 280.345
>>>>>>>>>> La1 : 36.405
>>>>>>>>>>
>>>>>>>>>> float lon_0 ( lon_0 )
>>>>>>>>>> long_name : longitude
>>>>>>>>>> grid_type : Latitude/Longitude
>>>>>>>>>> units : degrees_east
>>>>>>>>>> Dj : 0.01
>>>>>>>>>> Di : 0.009999676
>>>>>>>>>> Lo2 : 280.735
>>>>>>>>>> La2 : 36.695
>>>>>>>>>> Lo1 : 280.345
>>>>>>>>>> La1 : 36.405
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On Sun, May 28, 2017 at 12:48 AM, Barry Lynn <
>>>>>>>>>> barry.h.lynn at gmail.com> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi:
>>>>>>>>>>>
>>>>>>>>>>> When the data is skewed, it is often because the projection
>>>>>>>>>>> you're reading is not the one your plotting.
>>>>>>>>>>>
>>>>>>>>>>> Do you know the projection of the data you're using?
>>>>>>>>>>>
>>>>>>>>>>> If you print(a) where a is the file you specify from addfile,
>>>>>>>>>>> you can see the variables (and hopefully projection) indicated within.
>>>>>>>>>>>
>>>>>>>>>>> Barry
>>>>>>>>>>>
>>>>>>>>>>> On Sat, May 27, 2017 at 4:20 PM, Tim Melino <melino33 at gmail.com>
>>>>>>>>>>> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Hi,
>>>>>>>>>>>>
>>>>>>>>>>>> I am attempting to plot a subset of the MRMS data from NCEP.
>>>>>>>>>>>> The latitude and longitude I am plotting is :
>>>>>>>>>>>>
>>>>>>>>>>>> LAT: 36.5459089
>>>>>>>>>>>> LON: -79.4585847
>>>>>>>>>>>>
>>>>>>>>>>>> I attached an image of the output along with the grib file I am
>>>>>>>>>>>> using for input. I am noticing a couple things first the data appears to
>>>>>>>>>>>> be skewed, so I am not sure if it is an issue with the map projection.
>>>>>>>>>>>> Second the point which should be plotted in Virginia is actually plotted in
>>>>>>>>>>>> North Carolina, the only thing that I have been able to do if offset the
>>>>>>>>>>>> marker to compensate for the issue. Does anyone have any ideas on how to
>>>>>>>>>>>> fix these issues? Any ideas would be greatly appreciated!
>>>>>>>>>>>>
>>>>>>>>>>>> Thanks!
>>>>>>>>>>>> Tim
>>>>>>>>>>>>
>>>>>>>>>>>> Code:
>>>>>>>>>>>>
>>>>>>>>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
>>>>>>>>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>>>>>>>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>>>>>>>>>>
>>>>>>>>>>>> begin
>>>>>>>>>>>>
>>>>>>>>>>>> lati=LAT
>>>>>>>>>>>> loni=LON
>>>>>>>>>>>> date=DATE
>>>>>>>>>>>> address=ADDRESS
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> ; Convert the Lat/Lon to a Float
>>>>>>>>>>>> latc= stringtofloat(lati)
>>>>>>>>>>>> lonc = stringtofloat(loni)
>>>>>>>>>>>>
>>>>>>>>>>>> ; Set the Map Bounds
>>>>>>>>>>>> minlat = latc - .12
>>>>>>>>>>>> maxlat = latc + .12
>>>>>>>>>>>> minlon = lonc - .15
>>>>>>>>>>>> maxlon = lonc + .15
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> ;========================
>>>>>>>>>>>> ; get list of all files and open as "one big file"
>>>>>>>>>>>> ;========================
>>>>>>>>>>>> a = addfile( "./" + lati + "_" + loni + ".grib2","r")
>>>>>>>>>>>>
>>>>>>>>>>>> type = "png"
>>>>>>>>>>>> wks = gsn_open_wks(type,"./" + lati + "" + loni + ".png")
>>>>>>>>>>>>
>>>>>>>>>>>> ; Set some Basic Plot options
>>>>>>>>>>>> res = True
>>>>>>>>>>>> res at gsnDraw = False ; do not draw
>>>>>>>>>>>> the plot
>>>>>>>>>>>> res at gsnFrame = False ; do not advance
>>>>>>>>>>>> the frame
>>>>>>>>>>>> pltres = True
>>>>>>>>>>>> pltres at NoTitles = True
>>>>>>>>>>>> pltres at CommonTitle = True
>>>>>>>>>>>> pltres at PlotTitle = ""
>>>>>>>>>>>> pltres at PanelPlot = True
>>>>>>>>>>>> pltres at FramePlot = False
>>>>>>>>>>>>
>>>>>>>>>>>> ; Set the Colormap
>>>>>>>>>>>> gsn_define_colormap(wks,"hail")
>>>>>>>>>>>> lat = a->lat_0
>>>>>>>>>>>> lon = a->lon_0
>>>>>>>>>>>> water = a->VAR_209_3_34_P0_L102_GLL0(:,:)
>>>>>>>>>>>> hail = water *0.0393701;
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> ; Plotting options for Hail
>>>>>>>>>>>> res at gsnMaximize = True
>>>>>>>>>>>> res at mpGeophysicalLineThicknessF = 4.0
>>>>>>>>>>>> ;res at mpProjection = "LambertConformal"
>>>>>>>>>>>> res at mpGridLineThicknessF = 0.5
>>>>>>>>>>>> res at mpLimbLineThicknessF = 4.0
>>>>>>>>>>>> res at mpNationalLineThicknessF = 4.0
>>>>>>>>>>>> res at mpUSStateLineThicknessF = 4.0
>>>>>>>>>>>> res at mpProvincialLineThicknessF = 4.0
>>>>>>>>>>>> res at mpOutlineBoundarySets = "AllBoundaries"
>>>>>>>>>>>> res at mpDataBaseVersion = "Ncarg4_1"
>>>>>>>>>>>> res at mpDataSetName = "Earth..4"
>>>>>>>>>>>> res at mpFillOn = False ; fill continents
>>>>>>>>>>>> res at mpPerimOn = False ; no box around map
>>>>>>>>>>>> res at tiMainOn = True
>>>>>>>>>>>> res at tiMainFontColor = "Black"
>>>>>>>>>>>> res at tiMainString =""
>>>>>>>>>>>> res at tiXAxisString = ""
>>>>>>>>>>>> res at tiYAxisString = ""
>>>>>>>>>>>> res at gsnLeftString = "Date: " + date ; add
>>>>>>>>>>>> the gsn titles
>>>>>>>>>>>> res at gsnRightString = ""
>>>>>>>>>>>> res at gsnStringFontHeightF = 0.016
>>>>>>>>>>>> res at cnInfoLabelOrthogonalPosF = 0.07 ; offset second label
>>>>>>>>>>>> information
>>>>>>>>>>>> res at gsnContourLineThicknessesScale = 0.001
>>>>>>>>>>>> res at cnFillOn = True
>>>>>>>>>>>> res at lbTitleOn = True ;
>>>>>>>>>>>> remove field name from label bar
>>>>>>>>>>>>
>>>>>>>>>>>> ; Set the Map Bounds
>>>>>>>>>>>> res at sfYArray = lat
>>>>>>>>>>>> res at sfXArray = lon
>>>>>>>>>>>> res at mpLimitMode = "LatLon"
>>>>>>>>>>>> res at mpMinLatF = minlat
>>>>>>>>>>>> res at mpMaxLatF = maxlat
>>>>>>>>>>>> res at mpMinLonF = minlon
>>>>>>>>>>>> res at mpMaxLonF = maxlon
>>>>>>>>>>>> res at tmXBOn = False
>>>>>>>>>>>> res at tmYLOn = False
>>>>>>>>>>>> res at tmXTOn = False
>>>>>>>>>>>> res at tmYROn = False
>>>>>>>>>>>> res at gsnAddCyclic = False
>>>>>>>>>>>> res at tfDoNDCOverlay = False ;
>>>>>>>>>>>> do not transform data
>>>>>>>>>>>> res at cnLevelSelectionMode = "ExplicitLevels" ; set
>>>>>>>>>>>> explicit contour levels
>>>>>>>>>>>> res at cnLevels = (/ 0.5,0.75, 1.00, 1.25, 1.50,
>>>>>>>>>>>> 1.75, 2.00, 2.50, 3.00, 3.50, 4.00, 4.50, 5.00/) ; set levels
>>>>>>>>>>>> res at tiMainFontHeightF = 0.018 ; font
>>>>>>>>>>>> height of title
>>>>>>>>>>>> res at lbAutoManage = False ; we
>>>>>>>>>>>> control label bar
>>>>>>>>>>>> res at pmLabelBarDisplayMode = "Always" ; turns
>>>>>>>>>>>> on label bar
>>>>>>>>>>>> res at lbOrientation = "Horizontal" ; ncl
>>>>>>>>>>>> default is vertical
>>>>>>>>>>>> res at pmLabelBarSide = "Bottom" ;
>>>>>>>>>>>> default is right
>>>>>>>>>>>> res at lbLabelStride = 1 ; skip
>>>>>>>>>>>> every other label
>>>>>>>>>>>> res at pmLabelBarWidthF = 0.7 ;
>>>>>>>>>>>> default is shorter
>>>>>>>>>>>> res at pmLabelBarHeightF = 0.1 ;
>>>>>>>>>>>> default is taller
>>>>>>>>>>>> res at lbLabelFontHeightF = .012 ;
>>>>>>>>>>>> default is HUGE
>>>>>>>>>>>> res at lbPerimOn = False ;
>>>>>>>>>>>> default has box
>>>>>>>>>>>> res at lbTitleString = "Hail Size (Inches)"
>>>>>>>>>>>> res at lbTitleFontHeightF = 0.012
>>>>>>>>>>>> plot = gsn_csm_contour_map(wks,hail,res)
>>>>>>>>>>>>
>>>>>>>>>>>> ; Setup the point
>>>>>>>>>>>> mstring1 = "z"
>>>>>>>>>>>> fontnum1 = 35
>>>>>>>>>>>> xoffset1 = .4
>>>>>>>>>>>> yoffset1 = .7
>>>>>>>>>>>> ratio1 = 0.0
>>>>>>>>>>>> size1 = 1.0
>>>>>>>>>>>> angle1 = 0.0
>>>>>>>>>>>>
>>>>>>>>>>>> Tm_index = NhlNewMarker(wks, mstring1, fontnum1, xoffset1,
>>>>>>>>>>>> yoffset1, ratio1, size1, angle1)
>>>>>>>>>>>>
>>>>>>>>>>>> ; Plot the sites and label them
>>>>>>>>>>>> gres = True
>>>>>>>>>>>> gres at gsMarkerIndex = Tm_index ; 16
>>>>>>>>>>>> gres at gsMarkerSizeF = 0.015 ; .009
>>>>>>>>>>>> gres at gsMarkerColor = "black"
>>>>>>>>>>>> gres at txFont ="helvetica-bold"
>>>>>>>>>>>> gres at txFontHeightF = "25"
>>>>>>>>>>>>
>>>>>>>>>>>> site = (/address/)
>>>>>>>>>>>> lats = (/latc/)
>>>>>>>>>>>> lons = (/lonc/)
>>>>>>>>>>>>
>>>>>>>>>>>> dot = gsn_add_polymarker(wks,plot,lons,lats,gres)
>>>>>>>>>>>> tlats =lats + .03
>>>>>>>>>>>> tlons =lons
>>>>>>>>>>>> text = gsn_add_text(wks,plot,site,tlons,tlats ,gres)
>>>>>>>>>>>>
>>>>>>>>>>>> ; Draw the graphics
>>>>>>>>>>>> draw(plot)
>>>>>>>>>>>> frame(wks)
>>>>>>>>>>>>
>>>>>>>>>>>> end
>>>>>>>>>>>>
>>>>>>>>>>>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170601/bd135a8f/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sample_point.png
Type: image/png
Size: 95010 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170601/bd135a8f/attachment.png
More information about the ncl-talk
mailing list