[ncl-talk] Area average using a shapefile
Karin Meier-Fleischer
meier-fleischer at dkrz.de
Wed Apr 15 12:04:14 MDT 2020
Hi Lyndz,
in your case I think you can't use the shapefile_mask_data function
because you want to select an element of the variable Type. By the way
the shapefile is not correct because ID is to written properly.
I've attached the NCL script to compute the average of the polygon for
Type equal 2.000. For your data there is only one value inside the
shapefile polygon.
-Karin
Am 14.04.20 um 14:22 schrieb Lyndz:
> Hi Karin,
>
> Thank you for the response.
>
> The attached file is the input NetCDF file.
>
> Sincerely,
> Lyndz
>
> On Tue, Apr 14, 2020 at 8:11 PM Karin Meier-Fleischer
> <meier-fleischer at dkrz.de <mailto:meier-fleischer at dkrz.de>> wrote:
>
> Hi Lyndz,
>
> excuse me for the late reply. Without the data file
> ymonmean_temp_1979-2018.nc <http://ymonmean_temp_1979-2018.nc> it is
> not possible to help.
>
> Sorry,
> Karin
>
> Am 04.04.20 um 13:51 schrieb Lyndz:
> > Dear NCL-experts,
> >
> > *[Description of the Problem]*
> >
> > I would like to get the area average of temperature (netcdf file)
> based
> > on a shapefile mask.
> >
> > The attached file contains a shapefile of climate types.
> > I want to get the area average based on Climate type 2.
> > I want to do this per time step and save the output file to a
> text file.
> >
> >
> > *[What I have so far]*
> >
> > I have a script that plots the shapefile. However, I dont know
> how I can
> > get the area average of climate type 2 only.
> >
> >
> > *[Error]*
> > I get the following error using the attached script:
> >
> > (0)gsn_csm_contour_map_other: Fatal: the input data array must be
> 1D or 2D
> >
> > fatal:Illegal right-hand side type for assignment
> >
> > fatal:["Execute.c":8637]:Execute: Error occurred at or near line
> 65 in
> > file mask_rainfall.ncl
> >
> >
> >
> > I'll appreciate any help on how I can do this correctly in NCL.
> >
> > Sincerely,
> > Lyndz
>
--
Dipl. Geophys. Karin Meier-Fleischer
Visualization Group - NCL, CDO, Python
Application Support
Deutsches Klimarechenzentrum GmbH (DKRZ)
Bundesstrasse 45a - D20146 Hamburg - Germany
Phone: +49 (0)40 460094 126
Fax: +49 (0)40 460094 270
E-Mail: meier-fleischer at dkrz.de
URL: www.dkrz.de
Geschäftsführer: Prof. Dr. Thomas Ludwig
Sitz der Gesellschaft: Hamburg
Amtsgericht Hamburg HRB 39784
-------------- next part --------------
load "$HOME/NCL/shapefiles/shapefile_utils.ncl"
;load "./shapefile_utils.ncl"
begin
;-- shapefile name
shp_filename = "PHL_climatetype.shp"
;-- open file and read variable
f = addfile("ymonmean_temp_1979-2018.nc","r")
var = f->air(:,{850},:,:)
time = 0 ;-- set timestep to be plotted
;-- define sub-region, here for Germany
minlat = 8
maxlat = 20
minlon = 120
maxlon = 128
;-- open workstation
wks_type = "png"
wks_type at wkWidth = 1200
wks_type at wkHeight = 1200
wks = gsn_open_wks(wks_type,"plot_shapefiles_masking")
;-- 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 mpDataBaseVersion = "HighRes" ;-- use HighRes map
res at mpDataResolution = "Fine" ;-- we need a finer resolution
res at mpProjection = "Mercator" ;-- use Mercator projection
res at mpLimitMode = "Corners" ;-- map limit mode
res at mpLeftCornerLatF = minlat-1.0 ;-- min lat
res at mpRightCornerLatF = maxlat+1.0 ;-- max lat
res at mpLeftCornerLonF = minlon-1.0 ;-- min lon
res at mpRightCornerLonF = maxlon+1.0 ;-- max lon
res at mpCenterLatF = ((minlat+maxlat)*0.5)
res at mpCenterLonF = ((minlon+maxlon)*0.5)
res at mpGridLineColor = "grey40"
res at mpFillOn = True
res at mpOutlineOn = True
res at mpGeophysicalLineColor = "black"
res at mpOceanFillColor = (/ 0.824, 0.961, 1.0 /)
res at mpInlandWaterFillColor = (/ 0.824, 0.961, 1.0 /)
res at mpLandFillColor = (/ 0.7, 0.7, 0.7 /)
res at pmTickMarkDisplayMode = "Always"
res at tmXTOn = False
res at tiMainString = "Masked by shapefile"
res at cnFillMode = "RasterFill"
res at cnFillOn = True
res at cnLinesOn = False
res at cnLineLabelsOn = False
res at cnFillPalette = "BlueYellowRed"
res at cnLevelSelectionMode = "ManualLevels" ;-- manually set contour levels
res at cnMinLevelValF = 260.0 ;-- contour min lev
res at cnMaxLevelValF = 290.0 ;-- contour max lev
res at cnLevelSpacingF = 0.5 ;-- contour spacing;;
res at lbBoxMinorExtentF = 0.16 ;-- 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
;-- create plot
plot = gsn_csm_contour_map(wks,var(time,:,:),res)
;---------------------------
;-- read data from shapefile
;---------------------------
shpf = addfile(shp_filename,"r")
segments = shpf->segments
geometry = shpf->geometry
segsDims = dimsizes(segments)
geomDims = dimsizes(geometry)
lon = shpf->x
lat = shpf->y
nType = 2 ;-- Type = 2.000
;-- get global attributes
geom_segIndex = shpf at geom_segIndex
geom_numSegs = shpf at geom_numSegs
segs_xyzIndex = shpf at segs_xyzIndex
segs_numPnts = shpf at segs_numPnts
numFeatures = geomDims(0)
lat1d = f->lat
lon1d = f->lon
nlat = dimsizes(lat1d)
nlon = dimsizes(lon1d)
;-- add polylines from the shape files to the plot
polylines = new(segsDims(0),graphic)
;-- set resources for the shapefile polygon
plres = True
plres at gsLineColor = "white"
plres at gsLineThicknessF = 4.0
;-- create 1D array for values inside shapefile polygon
var1d = ndtooned(var(time,:,:))
var1d = var1d at _FillValue
;-- set start and end segment values for selected nType
startSegment = geometry(nType, geom_segIndex)
numSegments = geometry(nType, geom_numSegs)
segNum = 0
k = 0
do seg=startSegment, startSegment+numSegments-1
startPT = segments(seg, segs_xyzIndex)
endPT = startPT + segments(seg, segs_numPnts) - 1
;-- draw the segments
polylines(segNum) = gsn_add_polyline(wks, plot, lon(startPT:endPT),\
lat(startPT:endPT), plres)
;-- find the values inside the shapefile polygon
do n=0,nlat-1
do m=0,nlon-1
if(lat1d(n) .lt. min(lat(startPT:endPT)) .or. lat1d(n) .gt. max(lat(startPT:endPT)).or.\
lon1d(m) .lt. min(lon(startPT:endPT)) .or. lon1d(m) .gt. max(lon(startPT:endPT)))
continue
end if
if(gc_inout(lat1d(n), lon1d(m), lat(startPT:endPT), lon(startPT:endPT))) then
var1d(k) = var(time,n,m) ;-- values inside the shapefile polygon
end if
k = k + 1
end do
end do
segNum = segNum + 1
end do
count = dimsizes(ind(.not. ismissing(var1d))) ;-- number of values inside shapefile polygon
average = avg(var1d) ;-- compute the average value
print("")
print("--> values inside shapefile polygon = " + count)
print("")
print("--> average value = " + average)
print("")
;-- draw the plot
draw(plot)
frame(wks)
end
More information about the ncl-talk
mailing list