Hi everyone,
I have tried to put into practise what Rick said below, and I am not getting any errors, but I'm also getting any graphs. The script just keeps running as if it is doing something, but doesn't stop or output anything. If I remove the lines corresponding to the mask it works again. Can anyone see any problem with the script below? This is the text that appears on the command line when I run the script:
(0) ==================================================
(0) Shapefile: /home/cohara/IRL_adm/IRL_adm0.shp
(0) Areas of interest: the whole shapefile
(0) min_lat_chk: 51.4196
(0) max_lat_chk: 55.4504
(0) min_lon_chk: -10.6631
(0) max_lon_chk: -5.99361
(0) min_lat_data: 50
(0) max_lat_data: 56.99
(0) min_lon_data: -11
(0) max_lon_data: -5
(0) 47034 data values originally
(0) Will keep data values inside given shapefile areas
This is the script:
--------------------------------------------------------------------------------------------------------------------------
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/WRFUserARW.ncl"
load "./scripts/rainfall/for_map/shapefile_mask_data.ncl"
begin ;8
n=0
i=1981
print("Hi")
y=systemfunc("ls /home/cohara/RainfallData/for_map/timsum_*.nc") ;12
x=addfiles(y,"r")
do while(i.le.2013)
data=x[n]->var61(0,0,:,:)
shpfile="/home/cohara/IRL_adm/IRL_adm0.shp" ;16
opt=True
opt@debug=True
;opt@shape_names=(/"Ireland"/)
data_mask=shapefile_mask_data(data,shpfile,opt) ;20
plres=True
plres@gsEdgesOn=True
plres@gsLineColor="black" ;24
plres@gsEdgeColor="black"
plres@PanelPlot=True
resources=True
resources@gsnMaximize=True ;28
resources@gsnDraw=False
resources@gsnFrame=False
resources@tiMainString="Total Annual Rainfall " + sprinti("%0.4i",i)
resources@cnMonoLineColor=False ;Turn of the drawing of contour lines in 1 colour ;32
resources@tiMainFontColor="black"
resources@tiMainFont="times-roman"
resources@cnFillOn=True ;Turn on contour line fill
resources@cnMonoFillPattern=True ;Turn on single fill pattern (default) ;36
resources@cnMonoFillColor=False ;(default)
resources@cnMonoLineColor=False ;(default)
resources@cnLinesOn=False ;(default)
resources@tmYLOn=False ;Removes tick marks from left of y axis ;40
resources@tmYROn=False
resources@tmXTOn=False
resources@tmXBOn=False
resources@cnLineLabelsOn=False ;Removes contour line labels ;44
resources@cnInfoLabelOn=False ;Removes the informational label
resources@pmLabelBarDisplayMode="Always";Turns on the label bar
resources@lbLabelBarOn=True
resources@lbBoxLinesOn=False ;Removes lines from between the label bar colours ;48
resources@lbTitleOn=True ;Turns on label bar title
resources@lbTitleString="Rainfall (mm)"
resources@lbPerimOn=False ;Removes perimeter from label bar
resources@lbTitleFontColor="black" ;52
resources@lbTitleFont="times-roman"
resources@lbOrientation="Horizontal"
resources@lbLabelPosition="Bottom"
resources@lbTitlePosition="Top" ;56
;resources@cnExplicitLegendLabelsOn=True
;resources@lgLabelStrings=(/"4","5","6","7","8","9","10","11","12","13","14"/)
;resources@cnLineLabelStrings=resources@lbLabelStrings
resources@lbTitleDirection="Across" ;60
resources@lbTitleFontHeightF=0.02
resources@lbLabelFont="times-roman"
resources@lbLabelFontColor="black"
;resources@cnLevelSelectionMode="ManualLevels" ;Set the label bar levels manually ;64
;resources@cnLevelSpacingF=200
resources@cnLevelSelectionMode="ExplicitLevels"
resources@cnLevels=(/400,500,600,700,800,900,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000/)
resources@mpProjection="Mercator" ;Map Projection ;68
resources@gsnAddCyclic=False
resources@mpLimitMode="LatLon" ;The viewable portion of the map determined by latitude and longitude
resources@mpMinLatF=51
resources@mpMaxLatF=55.7 ;72
resources@mpMinLonF=-11
resources@mpMaxLonF=-5.2
resources@mpGridAndLimbOn=True
resources@mpDataSetName="Earth..4" ;The data set used (most up to date) ;76
resources@mpDataBaseVersion="MediumRes"
resources@mpDataResolution="Finest" ;Map render resolution
resources@mpOutlineOn=False
resources@mpFillOn=False ;80
;resources@vpWidthF=1
;resources@vpHeightF=0.7
wtype="png"
;wtype@wkWidth=2500 ;84
;wtype@wkHeight=2500
;sfiles=(/"/home/cohara/IRL_adm/IRL_adm0","/home/cohara/GBR_adm/GBR_adm0"/)+".shp"
xwks=gsn_open_wks(wtype,"Total_" + sprinti("%0.4i",i))
gsn_define_colormap(xwks,"precip4_11lev") ;88
plot=gsn_csm_contour_map(xwks,data,resources)
poly0=gsn_add_shapefile_polylines(xwks,plot,"/home/cohara/IRL_adm/IRL_adm1.shp",plres)
poly1=gsn_add_shapefile_polylines(xwks,plot,"/home/cohara/GBR_adm/GBR_adm2.shp",plres)
;poly2=gsn_add_shapefile_polylines(xwks,plot,"/home/cohara/IRL_adm/IRL_adm0.shp",plres) ;92
contour_mask=wrf_contour(x,xwks,data_mask,resources)
plot_mask=wrf_map_overlays(x,xwks,contour_mask,plres,resources)
draw(plot)
frame(xwks) ;96
n=n+1
i=i+1
end do
end ;100
Thanks for any help,
Ciara
> From: brownrig@ucar.edu
> Subject: Re: fatal:mask: dimension sizes of parameter 0 and parameter 1 do not match
> To: cohara1985@hotmail.com; ncl-talk@ucar.edu
> Date: Wed, 5 Mar 2014 08:37:41 -0700
>
> Hi Ciara,
>
> On line 17:
>
> data_mask=mask(data,shpfile,opt)
>
> I think you want to use function shapefile_mask_data(). This is an
> external function that can be obtained by looking at example #14 at
> the link below; that example is also one of performing masking in the
> fashion you are trying to do:
>
> http://www.ncl.ucar.edu/Applications/shapefiles.shtml
>
> I hope that helps…
>
>
> Rick
>
> On Wed, 5 Mar 2014 12:40:17 +0000
> Ciara O'Hara <cohara1985@hotmail.com> wrote:
> > Hi,
> >
> > I am running the script below, and getting the following error
> >message:
> >
> > ------------------------------------------------------------------------------------------------------------------------
> > fatal:mask: dimension sizes of parameter 0 and parameter 1 do not
> >match
> > fatal:["Execute.c":8128]:Execute: Error occurred at or near line 17
> >in file /home/cohara/scripts/rainfall/for_map/contour_map_script.ncl
> > ------------------------------------------------------------------------------------------------------------------------
> >
> > I am trying to use a shapefile to mask out the sea, so that my plot
> >has nice clean edges around the coast. I understand that it is a
> >latitude and longitude dimension size problem, but I'm not sure how
> >to fix it. I've run print_Var_Summary on one of my data variables
> >(they are all the same except for the date attached), here is the
> >output:
> >
> > ------------------------------------------------
> > ncl 7> printVarSummary(data)
> >
> > Variable: data
> > Type: float
> > Total Size: 188136 bytes
> > 47034 values
> > Number of Dimensions: 2
> > Dimensions and sizes: [lat | 234] x [lon | 201]
> > Coordinates:
> > lat: [ 50..56.99]
> > lon: [ -11.. -5]
> > Number Of Attributes: 4
> > height : 0
> > time : 19810930
> > table : 1
> > _FillValue : -9e+33
> > ------------------------------------------------
> >
> > Any suggestions would be appreciated,
> >
> > Best Regards,
> > Ciara
> >
> > ---------------------------------
> > script:
> > ---------------------------------
> >
> > 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"
> >
> > begin
> > n=0
> > i=1981
> > print("Hi")
> > y=systemfunc("ls
> >/home/cohara/RainfallData/for_map/timsum_*.nc")
> > x=addfiles(y,"r")
> > do while(i.le.2013)
> > data=x[n]->var61(0,0,:,:)
> > shpfile="/home/cohara/IRL_adm/IRL_adm0.shp"
> > opt=True
> > opt@return_mask=True
> > data_mask=mask(data,shpfile,opt)
> > data_land_mask=where(data_mask.eq.1,data,data@_FillValue)
> > copy_VarMeta(data,data_land_mask)
> > plres=True
> > plres@gsEdgesOn=True
> > plres@gsLineColor="black"
> > plres@gsEdgeColor="black"
> > resources=True
> > resources@gsnMaximize=True
> > resources@gsnDraw=False
> > resources@gsnFrame=False
> > resources@tiMainString="Total Annual Rainfall " +
> >sprinti("%0.4i",i)
> > resources@cnMonoLineColor=False ;Turn of the drawing
> >of contour lines in 1 colour
> > resources@tiMainFontColor="black"
> > resources@tiMainFont="times-roman"
> > resources@cnFillOn=True ;Turn on contour line
> >fill
> > resources@cnMonoFillPattern=True ;Turn on single fill
> >pattern (default)
> > resources@cnMonoFillColor=False ;(default)
> > resources@cnMonoLineColor=False ;(default)
> > resources@cnLinesOn=False ;(default)
> > resources@tmYLOn=False ;Removes tick marks
> >from left of y axis
> > resources@tmYROn=False
> > resources@tmXTOn=False
> > resources@tmXBOn=False
> > resources@cnLineLabelsOn=False ;Removes contour line
> >labels
> > resources@cnInfoLabelOn=False ;Removes the
> >informational label
> > resources@pmLabelBarDisplayMode="Always";Turns on the label
> >bar
> > resources@lbLabelBarOn=True
> > resources@lbBoxLinesOn=False ;Removes lines from
> >between the label bar colours
> > resources@lbTitleOn=True ;Turns on label bar
> >title
> > resources@lbTitleString="Rainfall (mm)"
> > resources@lbPerimOn=False ;Removes perimeter
> >from label bar
> > resources@lbTitleFontColor="black"
> > resources@lbTitleFont="times-roman"
> > resources@lbOrientation="Horizontal"
> > resources@lbLabelPosition="Bottom"
> > resources@lbTitlePosition="Top"
> > ;resources@cnExplicitLegendLabelsOn=True
> > ;resources@lgLabelStrings=(/"4","5","6","7","8","9","10","11","12","13","14"/)
> > ;resources@cnLineLabelStrings=resources@lbLabelStrings
> > resources@lbTitleDirection="Across"
> > resources@lbTitleFontHeightF=0.02
> > resources@lbLabelFont="times-roman"
> > resources@lbLabelFontColor="black"
> > ;resources@cnLevelSelectionMode="ManualLevels" ;Set the
> >label bar levels manually
> > ;resources@cnLevelSpacingF=0.2
> > resources@cnLevelSelectionMode="ExplicitLevels"
> > resources@cnLevels=(/0,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800,4000/)
> > resources@mpProjection="Mercator" ;Map Projection
> > resources@gsnAddCyclic=False
> > resources@mpLimitMode="LatLon" ;The viewable portion
> >of the map determined by latitude and longitude
> > resources@mpMinLatF=51
> > resources@mpMaxLatF=55.7
> > resources@mpMinLonF=-11
> > resources@mpMaxLonF=-5.2
> > resources@mpGridAndLimbOn=True
> > resources@mpDataSetName="Earth..4" ;The data set used
> >(most up to date)
> > resources@mpDataBaseVersion="MediumRes"
> > resources@mpDataResolution="Finest" ;Map render
> >resolution
> > resources@mpOutlineOn=False
> > resources@mpFillOn=False
> > ;resources@vpWidthF=1
> > ;resources@vpHeightF=0.7
> > wtype="png"
> > ;wtype@wkWidth=2500
> > ;wtype@wkHeight=2500
> > ;sfiles=(/"/home/cohara/IRL_adm/IRL_adm0","/home/cohara/GBR_adm/GBR_adm0"/)+".shp"
> > xwks=gsn_open_wks(wtype,"Total_" + sprinti("%0.4i",i))
> > gsn_define_colormap(xwks,"precip4_11lev")
> > plot=gsn_csm_contour_map(xwks,data_land_mask,resources)
> > poly0=gsn_add_shapefile_polylines(xwks,plot,"/home/cohara/IRL_adm/IRL_adm1.shp",plres)
> > poly1=gsn_add_shapefile_polylines(xwks,plot,"/home/cohara/GBR_adm/GBR_adm2.shp",plres)
> > poly2=gsn_add_shapefile_polylines(xwks,plot,"/home/cohara/IRL_adm/IRL_adm0.shp",plres)
> > draw(plot)
> > frame(xwks)
> > n=n+1
> > i=i+1
> > end do
> > end
> >
> >
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Mar 14 06:11:11 2014
This archive was generated by hypermail 2.1.8 : Fri Mar 14 2014 - 15:08:52 MDT