[ncl-talk] How to avoid Interpolation outside shapefile
Rick Brownrigg
brownrig at ucar.edu
Wed Jan 8 21:19:45 MST 2020
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*
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>
> --
>
> *Thanks & Regards,*
>
> *Ajay Bankar.*
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200108/7fda3513/attachment.html>
More information about the ncl-talk
mailing list