[ncl-talk] Plotting Lat/Lon Points on a map
Tim Melino
melino33 at gmail.com
Wed May 31 19:45:57 MDT 2017
Hi Dave and others,
Thanks for your help on this. I figured it was some issue with the
boundaries in NCL, I will try to build in the GADM shapefile plotting as
you suggested. It would be good if the NCL database could be updated to
something a little more accurate, we are using NCL to plot some high
resolution data that is very zoomed in. This is not the first time we have
seen this type of issue.
Tim
On Tue, May 30, 2017 at 9: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
>>>>>>>>>>>>
>>>>>>>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170531/978e4a8c/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/20170531/978e4a8c/attachment.png
More information about the ncl-talk
mailing list