[ncl-talk] Shapefile Problem with 6.2.1
Scott Capps
scott at allvertum.com
Tue Nov 11 16:00:43 MST 2014
Greetings,
I have just upgraded to version 6.2.1
(ncl_ncarg-6.2.1.Linux_Debian7.6_x86_64_nodap_gcc472.tar.gz) and am
happy with the many enhancements. I am especially interested in the
more efficient shapefile plotting abilities. However, I am getting
strange output when trying the new shapefile methods. Attached is a ps
file generated using the new method (Which has an issue). I can provide
the shapefile (Major U.S. highways) as well.
; THIS IS THE NEW METHOD WHICH PRODUCES STRANGE OUTPUT:
plres = True ; resources for polylines
plres at gsLineColor = "gray25"
plres at gsLineThicknessF = 3
f =
addfile(ncl_dir+"shapefiles/intrastate_shapefiles/intrstat.shp", "r")
; Open shapefile
plres at gsSegments = f->segments ;(:,0)
poly0 = gsn_add_polyline(wks, plot, f->x, f->y, plres)
;
; THIS OLD METHOD WORKS BUT IS VERY SLOW:
; f =
addfile(ncl_dir+"shapefiles/intrastate_shapefiles/intrstat.shp", "r")
; Open shapefile
; segments = f->segments
; geometry = f->geometry
; segsDims = dimsizes(segments)
; geomDims = dimsizes(geometry)
; geom_segIndex = f at geom_segIndex
; geom_numSegs = f at geom_numSegs
; segs_xyzIndex = f at segs_xyzIndex
; segs_numPnts = f at segs_numPnts
; numFeatures = geomDims(0)
; lines = new(segsDims(0),graphic) ; array to hold polylines
; plres = True ; resources for polylines
; plres at gsLineColor = "gray25"
; plres at gsLineThicknessF = 3
; lon = f->x
; lat = f->y
; segNum = 0 ; Counter for adding polylines
; do i=0, numFeatures-1
; startSegment = geometry(i, geom_segIndex)
; numSegments = geometry(i, geom_numSegs)
; do seg=startSegment, startSegment+numSegments-1
; startPT = segments(seg, segs_xyzIndex)
; endPT = startPT + segments(seg, segs_numPnts) - 1
; lines(segNum) = gsn_add_polyline(wks, plot,
lon(startPT:endPT), \
; lat(startPT:endPT), plres)
; segNum = segNum + 1
; end do
; end do
;
The full NCL script is attached. I can provide sample data upon request
(invariant_d01.nc, wrftest.nc and the shapefile).
Thank you,
Scott
-------------- next part --------------
; ***********************************************************
; USAGE:
; ncl 'fil_inv="./invariant_d01.nc"' 'dir_in="./"' 'fl_in="wrftest.nc"' 'dir_out="./"' 'fl_out="test_plot"' 'fl_typ="ps"' test.ncl
; ***********************************************************
; ***********************************************************
; ***********************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCL_FXN_PATH/time_conversion_module.ncl"
begin
;
; search for existing model output at this time
file_in = addfile(dir_in+fl_in,"r")
; *************************************
; GET FIXED OUTPUT
; *************************************
init_time = stringtochar(file_in at SIMULATION_START_DATE)
dx = file_in at DX
cen_lon = file_in at CEN_LON
cen_lat = file_in at CEN_LAT
true_lat1 = file_in at TRUELAT1
true_lat2 = file_in at TRUELAT2
stand_lon = file_in at STAND_LON
; INVARIANT FILE
file1 = addfile(fil_inv,"r")
lat2d = file1->XLAT(0,:,:)
lon2d = file1->XLONG(0,:,:)
lat2d at units= "degrees_north"
lon2d at units= "degrees_east"
dimll = dimsizes(lat2d)
nlat = dimll(0)
mlon = dimll(1)
; *************************************
; END: GET FIXED OUTPUT
; *************************************
;
; 2m Temperature
left_title = "2m Temp (~S~o~N~F)/ Wind Barbs (Knots)"
color_map = "BlAqGrYeOrReVi200"
cnspc = 5
lblstrd = 1
lbltitle = "(~S~o~N~F)"
plot_var = file_in->T2
plot_var = (plot_var-273.15)*1.8+32.
plot_var at units = "F"
plot_var at lat2d = lat2d
plot_var at lon2d = lon2d
; GET MAX AND MIN FOR PLOTTING
max_tpt = 115.
min_tpt = -10.
; WSPD Components 10m (m/s) Earth Relative
u10m = file_in->U10
v10m = file_in->V10
; Convert to knots for wind barbs
u10m = u10m * 1.943846
v10m = v10m * 1.943846
delete(file_in)
;
u10m at units = "knots"
v10m at units = "knots"
u10m at lat2d = lat2d
u10m at lon2d = lon2d
v10m at lat2d = lat2d
v10m at lon2d = lon2d
;
nstart = 0
nstride = 1
;
wks = gsn_open_wks("ps","test.ps")
gsn_define_colormap(wks,color_map)
gsn_define_colormap(wks,color_map)
gsn_reverse_colormap(wks)
;
plot = new(1,graphic)
res = True
res at pmTickMarkDisplayMode = "NoCreate"
res at gsnMaximize = True
res at gsnDraw = False ; don't draw
res at gsnFrame = False ; don't advance frame
res at mpFillOn = False
;
res at mpProjection = "LambertConformal"
res at mpLambertParallel1F = get_res_value_keep(res, "mpLambertParallel1F",file1 at TRUELAT1)
res at mpLambertParallel2F = get_res_value_keep(res, "mpLambertParallel2F",file1 at TRUELAT2)
res at mpLambertMeridianF = get_res_value_keep(res, "mpLambertMeridianF",file1 at STAND_LON)
;
; Full domain
; 867 wide x 872 tall
res at mpLimitMode = "Corners"
res at mpLeftCornerLatF = lat2d(0,0)
res at mpLeftCornerLonF = lon2d(0,0)
res at mpRightCornerLatF = lat2d(nlat-1,mlon-1)
res at mpRightCornerLonF = lon2d(nlat-1,mlon-1)
vcmindist = 0.02
;
res at pmLabelBarOrthogonalPosF = .01 ; move whole thing down
; END: Regional settings
; ;
res at mpOutlineBoundarySets = "GeophysicalAndUSStates"; turn on states
res at mpDataSetName = "Earth..2"
res at mpDataBaseVersion = "MediumRes"
res at mpInlandWaterFillPattern = 2
res at mpUSStateLineThicknessF = 2
res at mpNationalLineThicknessF = 2
res at mpGeophysicalLineThicknessF = 2
res at gsnAddCyclic = False
;
res at cnInfoLabelOn = False
res at cnFillOn = True
res at cnLinesOn = False
res at cnLineLabelsOn = False
res at gsnSpreadColors = True
res at lbLabelBarOn = True
res at lbTitleOn = False
res at lbTitleFontHeightF = .014
res at lbOrientation = "vertical"
res at lbTitlePosition = "Right"
res at lbTitleDirection = "Across"
res at cnMinLevelValF = min_tpt
res at cnMaxLevelValF = max_tpt
res at cnLevelSpacingF = cnspc
res at lbLabelStride = lblstrd
res at cnFillDrawOrder = "Draw" ; draw contours first
;res at pmLabelBarOrthogonalPosF = .01 ; move whole thing down
;
res at cnLevelSelectionMode = "ManualLevels"
res at gsnLeftString = ""
res at gsnCenterString = ""
res at gsnLeftStringFontHeightF = 0.012
res at gsnRightStringFontHeightF = 0.010
;
;VECTOR SETTINGS
res at vcRefAnnoOn = False
res at vcGlyphStyle ="WindBarb"
res at vcRefAnnoOn = False
res at vcRefLengthF = 0.016 ; define length of vec ref
res at vcMinDistanceF = vcmindist ;0.01;0.02
res at vcVectorDrawOrder = "PostDraw" ; draw vectors las
res at vcRefAnnoArrowLineColor = "black" ; change ref vector color
res at vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref
res at gsnScalarContour = True
;
plot = gsn_csm_vector_scalar_map(wks,u10m(0,:,:),v10m(0,:,:),plot_var(0,:,:),res)
;
; THIS IS THE NEW METHOD WHICH PRODUCES STRANGE OUTPUT:
plres = True ; resources for polylines
plres at gsLineColor = "gray25"
plres at gsLineThicknessF = 3
f = addfile(ncl_dir+"shapefiles/intrastate_shapefiles/intrstat.shp", "r") ; Open shapefile
plres at gsSegments = f->segments ;(:,0)
poly0 = gsn_add_polyline(wks, plot, f->x, f->y, plres)
;
; THIS OLD METHOD WORKS BUT IS VERY SLOW:
; f = addfile(ncl_dir+"shapefiles/intrastate_shapefiles/intrstat.shp", "r") ; Open shapefile
; segments = f->segments
; geometry = f->geometry
; segsDims = dimsizes(segments)
; geomDims = dimsizes(geometry)
; geom_segIndex = f at geom_segIndex
; geom_numSegs = f at geom_numSegs
; segs_xyzIndex = f at segs_xyzIndex
; segs_numPnts = f at segs_numPnts
; numFeatures = geomDims(0)
; lines = new(segsDims(0),graphic) ; array to hold polylines
; plres = True ; resources for polylines
; plres at gsLineColor = "gray25"
; plres at gsLineThicknessF = 3
; lon = f->x
; lat = f->y
; segNum = 0 ; Counter for adding polylines
; do i=0, numFeatures-1
; startSegment = geometry(i, geom_segIndex)
; numSegments = geometry(i, geom_numSegs)
; do seg=startSegment, startSegment+numSegments-1
; startPT = segments(seg, segs_xyzIndex)
; endPT = startPT + segments(seg, segs_numPnts) - 1
; lines(segNum) = gsn_add_polyline(wks, plot, lon(startPT:endPT), \
; lat(startPT:endPT), plres)
; segNum = segNum + 1
; end do
; end do
;
draw(plot)
frame(wks)
delete(plot_var)
delete(lat2d)
delete(lon2d)
end
exit
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.png
Type: image/png
Size: 361871 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20141111/99787f8c/attachment-0001.png
More information about the ncl-talk
mailing list