[ncl-talk] Creating Canada Mask...too slow ):

[Rick Brownrigg]

Hi Michael,

As noted, there are 24519 individual loops/polygons that make up the border
with all its islands, big and small.
One thing you might try is to skip testing against the ones that are
composed of small numbers of points.

If you open that file in NCL interactively, and type:

f = addfile(.....)
print(f->segments(:,1))    ; number of x/y coord pairs in each polygon

you'll get a sense for the distribution. There are 198 with over 1000
points; depending upon your data and its resolution, I speculate you might
still get a good result skipping over anything with less than 1000;
obviously this is something to experiment with.

Hope that helps...
Rick


On Wed, Jun 5, 2019 at 9:21 AM Michael Notaro <mnotaro at wisc.edu> wrote:

> The script below reads in a Canadian GADM shapefile and
> creates a 0.1-degree gridded mask where canada=1 for grid cells within
> Canada
> and canada=0 for grid cells outside of Canada.  It works but probably will
> take
> an estimated 15 hours to run.  Can anyone recommend a way to speed it up?
> Thanks, Michael
>
>
> 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 "./shapefile_utils.ncl"
> begin
>
> a=addfile("create_gridded_snwd.nc","r")
> lat=a->lat
> lon=a->lon
>
> b=addfile("gadm36_CAN_shp/gadm36_CAN_0.shp","r") ; num_segments=24519
> shape_lon=b->x
> shape_lat=b->y
> geometry=b->geometry
> segments=b->segments
> segsDims = dimsizes(segments)
> geomDims = dimsizes(geometry)
> geom_segIndex = b at geom_segIndex
> geom_numSegs  = b at geom_numSegs
> segs_xyzIndex = b at segs_xyzIndex
> segs_numPnts  = b at segs_numPnts
> numFeatures = geomDims(0)
>
> canada=new((/201,301/),integer)
> canada=0
>
> segNum = 0
> 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
>     lons=shape_lon(startPT:endPT)
>     lats=shape_lat(startPT:endPT)
>     do ilat=60,200
>       do ilon=0,300
>         inout=gc_inout(lat(ilat),lon(ilon),lats,lons)
>         if (inout) then
>           canada(ilat,ilon)=1
>         end if
>       end do
>     end do
>     segNum = segNum + 1
>     delete(lons)
>     delete(lats)
>   end do
> end do
>
> canada!0="lat"
> canada!1="lon"
> canada&lat=lat
> canada&lon=lon
>
> system("rm create_canada_mask.nc")
> out=addfile("create_canada_mask.nc","c")
> out->canada=canada
> out->lat=lat
> out->lon=lon
>
> end
>
>
>
>
> Michael Notaro
> Associate Director
> Nelson Institute Center for Climatic Research
> University of Wisconsin-Madison
> Phone: (608) 261-1503
> Email: mnotaro at wisc.edu
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190605/26ab3033/attachment.html>