<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"><!-- P {margin-top:0;margin-bottom:0;} --></style>
</head>
<body dir="ltr">
<div id="divtagdefaultwrapper" style="font-size:12pt;color:#000000;font-family:Calibri,Helvetica,sans-serif;" dir="ltr">
<p style="margin-top:0;margin-bottom:0"></p>
<div>The script below reads in a Canadian GADM shapefile and</div>
<div>creates a 0.1-degree gridded mask where canada=1 for grid cells within Canada</div>
<div>and canada=0 for grid cells outside of Canada. It works but probably will take</div>
<div>an estimated 15 hours to run. Can anyone recommend a way to speed it up?</div>
<div>Thanks, Michael</div>
<div><br>
</div>
<div><br>
</div>
<div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"</div>
<div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"</div>
<div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"</div>
<div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"</div>
<div>load "./shapefile_utils.ncl"</div>
<div>begin</div>
<div><br>
</div>
<div>a=addfile("create_gridded_snwd.nc","r")</div>
<div>lat=a->lat</div>
<div>lon=a->lon</div>
<div><br>
</div>
<div>b=addfile("gadm36_CAN_shp/gadm36_CAN_0.shp","r") ; num_segments=24519</div>
<div>shape_lon=b->x</div>
<div>shape_lat=b->y</div>
<div>geometry=b->geometry</div>
<div>segments=b->segments</div>
<div>segsDims = dimsizes(segments)</div>
<div>geomDims = dimsizes(geometry)</div>
<div>geom_segIndex = b@geom_segIndex</div>
<div>geom_numSegs = b@geom_numSegs</div>
<div>segs_xyzIndex = b@segs_xyzIndex</div>
<div>segs_numPnts = b@segs_numPnts</div>
<div>numFeatures = geomDims(0)</div>
<div><br>
</div>
<div>canada=new((/201,301/),integer)</div>
<div>canada=0</div>
<div><br>
</div>
<div>segNum = 0</div>
<div>do i=0, numFeatures-1</div>
<div> startSegment = geometry(i, geom_segIndex)</div>
<div> numSegments = geometry(i, geom_numSegs)</div>
<div> do seg=startSegment, startSegment+numSegments-1</div>
<div> startPT = segments(seg, segs_xyzIndex)</div>
<div> endPT = startPT + segments(seg, segs_numPnts) - 1</div>
<div> lons=shape_lon(startPT:endPT)</div>
<div> lats=shape_lat(startPT:endPT)</div>
<div>
<div> do ilat=60,200</div>
<div> do ilon=0,300</div>
<div> inout=gc_inout(lat(ilat),lon(ilon),lats,lons)</div>
<div> if (inout) then</div>
<div> canada(ilat,ilon)=1</div>
<div> end if</div>
<div> end do</div>
<div> end do</div>
<div> segNum = segNum + 1</div>
<div> delete(lons)</div>
<div> delete(lats)</div>
<div> end do</div>
<div>end do</div>
<div><br>
</div>
<div>canada!0="lat"</div>
<div>canada!1="lon"</div>
<div>canada&lat=lat</div>
<div>canada&lon=lon</div>
<div><br>
</div>
<div>system("rm create_canada_mask.nc")</div>
<div>out=addfile("create_canada_mask.nc","c")</div>
<div>out->canada=canada</div>
<div>out->lat=lat</div>
<div>out->lon=lon</div>
<div><br>
</div>
<div>end</div>
<div><br>
</div>
<br>
</div>
<br>
<p></p>
<p style="margin-top:0;margin-bottom:0"><br>
</p>
<div id="Signature">
<div id="divtagdefaultwrapper" style="font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255); font-family: Calibri, Arial, Helvetica, sans-serif, Helvetica, EmojiFont, "Apple Color Emoji", "Segoe UI Emoji", NotoColorEmoji, "Segoe UI Symbol", "Android Emoji", EmojiSymbols;">
Michael Notaro<br>
<div>Associate Director<br>
</div>
<div>Nelson Institute Center for Climatic Research<br>
</div>
<div>University of Wisconsin-Madison<br>
</div>
<div>Phone: (608) 261-1503<br>
</div>
<div>Email: mnotaro@wisc.edu<br>
</div>
</div>
</div>
</div>
</body>
</html>