;;Usage: plot permafrost extent and label the research area
;;Author: Heyl 20150816


wks=gsn_open_wks("eps","permafrost")

res=True
res@gsnFrame=False
res@gsnDraw=False

res@mpProjection="LambertEqualArea"

;res@mpLimitMode="Corners"
;res@mpLeftCornerLatF=80
;res@mpLeftCornerLonF=50
;res@mpRightCornerLatF=20
;res@mpRightCornerLonF=150
res@mpLimitMode="LatLon"
res@mpMinLatF=20
res@mpMaxLatF=80
res@mpMinLonF=30
res@mpMaxLonF=150

res@pmTickMarkDisplayMode="Always"
res@tiMainString="Permafrost Extent"
plot=gsn_csm_map(wks,res)
print("ok")
;;;;;;;;;;;;;;;;;;
;;section to add polygons to map
;;;;;;;;;;;;;;;;;;
pathin="../original/"
f=addfile(pathin+"ggd600_permext_russia.shp","r")
;
;read data of shapefile
segments=f->segments
geometry=f->geometry
segsDims=dimsizes(segments)
geomDims=dimsizes(geometry)

geom_segIndex=f@geom_segIndex
geom_numSegs=f@geom_numSegs
segs_xyzIndex=f@segs_xyzIndex
segs_numPnts=f@segs_numPnts

lines=new(segsDims(0),graphic)
numFeatures=geomDims(0)

extent=f->EXTENT_OF_
lon=f->x
lat=f->y

plres=True
plres@gsEdgesOn=True
plres@gsEdgeColor="black"

colors=(/"blue","burlywood","azure","aquamarine","darkgreen","blue"/)
category=(/"Baikal","ground ice of seasonal perm.ro","<50","50-90",">90"/)
segNum=0
print("ok2")
do i=0,numFeatures-1
if (extent(i) .eq. category(0)) then
plres@gsFillColor=colors(0)
end if
if (extent(i) .eq. category(1)) then
plres@gsFillColor=colors(1)
end if
if (extent(i) .eq. category(2)) then
plres@gsFillColor=colors(2)
end if
if (extent(i) .eq. category(3)) then
plres@gsFillColor=colors(3)
end if
if (extent(i) .eq. category(4)) then
plres@gsFillColor=colors(4)
end if

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_polygon(wks,plot,lon(startPT:endPT),lat(startPT:endPT),plres)
segNum=segNum+1
end do
end do
print("ok3")
labels=(/"Water","Seasonal","<50 sporadic","50-90",">90 continuous"/)

lbres=True
lbres@vpWidthF=0.7
lbres@vpHeightF=0.05
lbres@lbPerimOn=False
lbres@lbOrientation="Horizontal"
lbres@lbLabelAlignment="BoxCenters"
lbres@lbFillColors=colors
lbres@lbMonoFillPattern=True
lbres@lbLabelFontHeightF=0.012
lbres@lbTitlePosition="Bottom"
lbres@lbTitleFontHeightF=0.01
gsn_labelbar_ndc(wks,5,labels,0.23,0.15,lbres)
;maximize_output(wks,False)

draw(plot)
frame(wks)