[ncl-talk] How to avoid Interpolation outside shapefile
Dave Allured - NOAA Affiliate
dave.allured at noaa.gov
Wed Jan 8 23:14:58 MST 2020
NCL can do both types of masking: data masking, and graphics clipping.
The best example I have seen is shapefile_21 on the Examples page. Look
closely at figures 2 and 3 to appreciate the differences. IMO, for
publication quality plots, graphics clipping is much superior to data
masking. Caveats for graphics clipping using a shapefile:
* The shapefile must be polygon type, not polyline.
* The shapefile must include surrounding areas on all sides. This is
because NCL is missing one or two key functions to support direct clipping
by shapefile. It uses a trick instead. The trick is to select all of the
areas that you do NOT want, and color them white (or background). This
lets the area of interest shine through the remaining holes. This is all
shown in the code for shapefile_21.
* If you have a full data field, do not use data masking underneath
graphics clipping. Allow the full data field to interpolate naturally.
Data masking can leave visible jagged edges and interpolation distortion
around the edges.
There are several ways to obtain a suitable polygon shapefile for graphics
clipping. Extracting from a shapefile database is easiest, if you can find
the right database. A single region shapefile may be inverted into a
suitable "negative" polygon shapefile using good graphics software, or by
manual conversion, which is tedious.
On Wed, Jan 8, 2020 at 9:28 PM Rick Brownrigg via ncl-talk <
ncl-talk at ucar.edu> wrote:
> Ajay,
>
> Thanks for sending the files. I experimented with this for a good while
> and learned something I was not aware of. First I'd point out that when
> working with irregularly spaced point data, NCL must employ a triangulation
> of the data, and it is not uncommon to see strange-ish linear contour lines
> along the outer edges of a resultant contour plot, particularly if the
> point data has a concave shape, which your's does.
>
> So what I realized is that shapefile_mask_data() works in *data* space,
> not in *graphics* space (i.e., contour lines). All of the station points in
> your dataset are obviously inside the boundaries of the shapefile, so it
> returns a mask that is all 1's (inside). Then the contouring algorithm runs
> against that result, which in effect is the entire original concave
> dataset. The inevitable "strange-ish" contour lines form, but they are in
> graphics space, and hence, they are not clipped against the shapefile.
>
> So...I'm not sure what to tell you. One possibility would be to
> interpolate a regular grid from the station data (I plotted the individual
> stations and they seem to be quite dense, so that the resultant grid might
> make a good representation). That said, I don't know offhand what
> facilities are available for doing that -- perhaps that's a separate post
> to ncl-talk.
>
> Wish I had a better answer, but I hope that helps explain what it going on.
>
> Rick
>
>
> On Wed, Jan 8, 2020 at 8:03 PM Ajay Bankar <ajaybankar123 at gmail.com>
> wrote:
>
>> Hi Rick,
>> The same shapefile works fine with GPM data and wrf output
>> data but only having issue with station data. Please find attached
>> shapefiles. I have attached one plot in which I have used same shapefile
>> for wrf model output and GPM data. Here plots (a,b,c) represents wrf plots
>> and (d) represent GPM data plot and (e) represent station data plot.
>>
>> This is the output of ncl_filedump District_Boundary.shp
>>
>> Variable: f
>> Type: file
>> filename: District_Boundary
>> path: District_Boundary.shp
>> file global attributes:
>> layer_name : District_Boundary
>> geometry_type : polygon
>> geom_segIndex : 0
>> geom_numSegs : 1
>> segs_xyzIndex : 0
>> segs_numPnts : 1
>> dimensions:
>> geometry = 2
>> segments = 2
>> num_features = 30 // unlimited
>> num_segments = 33
>> num_points = 164066
>> variables:
>> integer geometry ( num_features, geometry )
>>
>> integer segments ( num_segments, segments )
>>
>> double x ( num_points )
>>
>> double y ( num_points )
>>
>> double PERIMETER ( num_features )
>>
>> double AREA ( num_features )
>>
>> int64 ID ( num_features )
>>
>> string DIST_NAME ( num_features )
>>
>> double SQ_KMS ( num_features )
>>
>>
>> On Thu, Jan 9, 2020 at 4:50 AM Rick Brownrigg <brownrig at ucar.edu> wrote:
>>
>>> Can you either i) provide the shapefile (and I'll point out that means
>>> the .shp, .shx, .dbf files), or perhaps show us the results of
>>>
>>> ncl_filedump District_Boundary.shp
>>>
>>> Rick
>>>
>>>
>>> On Wed, Jan 8, 2020 at 7:24 AM Ajay Bankar via ncl-talk <
>>> ncl-talk at ucar.edu> wrote:
>>>
>>>> Dear NCL Users,
>>>> I'm trying to make spatial plot of rainfall
>>>> from station data available in csv file. I tried different options in
>>>> resources but data is not masking with shapefile.
>>>> I have attached plot and data file for reference.
>>>>
>>>> Thanks for any help!
>>>>
>>>> Here is my script
>>>>
>>>> load "/shp.ncl"
>>>>
>>>> begin
>>>> ;-- shapefile name
>>>> shp_filename =
>>>> "/scratch/Administrative_Boundary_Headquarters/District_Boundary.shp"
>>>>
>>>> infile = "station_data.csv"
>>>> lines = asciiread(infile, -1, "string")
>>>>
>>>> rain = tofloat( str_get_field(lines(1:), 7, ","))
>>>> rain at lat1d = tofloat( str_get_field(lines(1:), 5, ","))
>>>> rain at lon1d = tofloat( str_get_field(lines(1:), 6, ","))
>>>> ;print(rain)
>>>>
>>>> ;-- define sub-region, here for Karnataka
>>>> minlat = 11
>>>> maxlat = 19
>>>> minlon = 73.5
>>>> maxlon = 79.1
>>>>
>>>> ;-- open workstation
>>>> wks_type = "pdf"
>>>> wks = gsn_open_wks(wks_type, "station_plot")
>>>> ;-- resource settings
>>>> res = True
>>>> res at gsnDraw = False ;--
>>>> don't draw plot yet
>>>> res at gsnFrame = False ;--
>>>> don't advance frame yet
>>>> res at gsnMaximize = True ;--
>>>> maximize plot in frame
>>>>
>>>> ; res at tfDoNDCOverlay = True
>>>>
>>>> cmap := read_colormap_file("BlAqGrYeOrReVi200")
>>>> cmap(0,:) = (/0,0,0,0/) ; make first color fully transparent
>>>> res at cnFillOn = True
>>>> res at cnLinesOn = False
>>>> res at cnLineLabelsOn = False
>>>> res at cnFillPalette = cmap
>>>> ;"BlueYellowRed"
>>>> res at cnLevelSelectionMode = "ExplicitLevels" ;--
>>>> manually set contour levels
>>>> res at cnLevels = (/5,20,35,50,65,80,95,110,125,140,155/)
>>>> res at gsnAddCyclic = False
>>>>
>>>> res at mpGridLineColor = 0 ;"grey40"
>>>> res at mpFillOn = False
>>>> res at mpOutlineOn = False
>>>> res at mpGeophysicalLineColor = "black"
>>>> res at mpLimitMode = "LatLon"
>>>> res at mpMinLatF = minlat
>>>> res at mpMaxLatF = maxlat
>>>> res at mpMinLonF = minlon
>>>> res at mpMaxLonF = maxlon
>>>> ;res at mpAreaMaskingOn = True
>>>>
>>>> res at lbBoxMinorExtentF = 0.1 ;--
>>>> decrease height of labelbar boxes
>>>> res at pmLabelBarOrthogonalPosF = 0.08 ;-- move
>>>> labelbar to the left side of plot
>>>> res at lbBoxMinorExtentF = 0.2 ;--
>>>> decrease height of labelbar boxes
>>>>
>>>> res at pmTickMarkDisplayMode = "Always"
>>>> res at tmXTOn = False
>>>> res at tiMainString = "STATION DATA PLOT"
>>>>
>>>> ;-- create plot
>>>> plot = gsn_csm_contour_map(wks,rain,res)
>>>> var_mask = shapefile_mask_data(rain,shp_filename,True)
>>>>
>>>> ;-- create contours of masked data
>>>> plot_mask = gsn_csm_contour_map(wks,var_mask,res)
>>>>
>>>> ;-- add plot_mask to plot
>>>> lnres = True
>>>> lnres at gsLineColor = "black"
>>>> lnres at gsLineThicknessF = 1.0
>>>>
>>>> plot at lines = gsn_add_shapefile_polylines(wks, plot_mask,
>>>> shp_filename, lnres)
>>>>
>>>> ;-- draw the plot and advance the frame
>>>> draw(plot_mask)
>>>> frame(wks)
>>>> end
>>>> --
>>>>
>>>> *Thanks & Regards,*
>>>> *Ajay*
>>>>
>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200108/40e9e1b4/attachment.html>
More information about the ncl-talk
mailing list