<div dir="ltr"><div><div><div><div><div><div><div><div><div><div><div>Hi Stavros,<br><br></div>I'm replying back to the list as the technique I'm suggesting may be useful to others.<br><br></div>As I understand your original problem was that masking was taking a long time, and you need to apply that make to a large number of datasets. I found that with the shapefile you sent, it was taking hours, not minutes. A simpler, more generalized shapefile would go a long ways towards alleviating the problem. But, given the shapefile at hand...<br><br></div>My thought was to create the mask once, and write it to a NetCDF file. Then when processing the individual files, load that mask and apply it against the data. The attached script shows most of the steps involved (it does not write to NetCDF, but that's simple enough -- examples at <a href="http://ncl.ucar.edu/Applications/o-netcdf.shtml" target="_blank">http://ncl.ucar.edu/<wbr>Applications/o-netcdf.shtml</a>).<br><br></div>- reads the lat/lon and snow data. Gets the x/y from the shapefile.<br></div>- creates a boolean "mask" variable the same dimensions as the snow data<br></div>- Loops over each of the individual polygons from the shapefile, and calls gc_inout() to find what part of the data grid overlaps individual polygons. gc_inout() returns a logical array denoting which grid cells are in/out. Those results are logically OR'd into the mask. On the premise that the smaller islands are not of concern relative to the grid resolution, the in/out test is skipped for polygons that have fewer vertices than a threshold (arbitrarily, 10,000 here -- that dragged in only the mainland and 3 other larger islands).  Even so, this loop/mask-generation phase takes 2-3hrs with the given shapefile.<br></div>- Use the mask and the where() function to mask out the data:<br><br>    snow = where(snow_mask.eq.True, snow, snow@_FillValue)<br><br></div>This assumes your data don't change spatial extent from file to file.  If this approach is useful, you'd want to break the script into two pieces: i) to generate and save the mask, ii) read your data, apply the mask, and generate the plots.<br><br></div>BTW: as a side comment, I may have thrown noise into the previous discussion regarding map projections. Your data and shapefile were both in geographic coordinates.<br><br></div>Hope this helps...<br></div>Rick<br><div><div><div><div><div><div><div><div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 1, 2017 at 12:27 PM, Stavros Dafis <span dir="ltr"><<a href="mailto:sdafis@cc.uoi.gr" target="_blank">sdafis@cc.uoi.gr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hello Rick,<br>
<br>
Thank you for answering my question about masking data in NCL. I am sorry for<br>
not replying earlier but I experienced some problems with my machines and I<br>
wanted to install Proj.4 and try your solution.<br>
<br>
I have not succeeded in converting the projection of my data to the projection<br>
of the shapefile. Maybe I am using Proj.4 in a wrong way and I found no<br>
particular documentation about the function "transform_coordinate".<br>
<br>
My data have the following projection:<br>
<br>
"Polar stereographic ellipsoidal projection with WGS-84 ellipsoid"<br>
<br>
I have uploaded an example of the data in netCDF, my NCL script and the<br>
shapefile of Greece:<br>
<a href="https://www.dropbox.com/sh/rfp5yxk1abfmmgf/AACGM3BfP26G9VgMUQeuD1kDa?dl=0" rel="noreferrer" target="_blank">https://www.dropbox.com/sh/rfp<wbr>5yxk1abfmmgf/AACGM3BfP26G9VgMU<wbr>QeuD1kDa?dl=0</a><br>
<br>
Any help would be very much appreciated.<br>
<br>
--<br>
Stavros NTAFIS (DAFIS)<br>
------------------------------<wbr>-----------------<br>
Physicist - Meteorologist, M.Sc.<br>
Ph.D. Candidate, Polytechnic School of Paris, France<br>
Laboratory of Dynamic Meteorology (LMD)<br>
Research Associate, National Observatory of Athens (NOA/IERSD)<br>
Tel. : <a href="tel:%28%2B33%291%2069%2033%2051%2064" value="+33169335164" target="_blank">(+33)1 69 33 51 64</a><br>
<a href="tel:%28%2B33%299%2081%2094%2022%2012" value="+33981942212" target="_blank">(+33)9 81 94 22 12</a><br>
Mobile: <a href="tel:%28%2B33%29%20066%20030%200147" value="+33660300147" target="_blank">(+33) 066 030 0147</a><br>
<a href="tel:%28%2B30%29%20697%2004%2020%20242" value="+306970420242" target="_blank">(+30) 697 04 20 242</a><br>
---------------<br>
Weather charts:<br>
<a href="http://www.meteo.gr" rel="noreferrer" target="_blank">http://www.meteo.gr</a><br>
<a href="http://www.meteo.gr/meteomaps/" rel="noreferrer" target="_blank">http://www.meteo.gr/meteomaps/</a><br>
<br>
<br>
Quoting Rick Brownrigg <<a href="mailto:brownrig@ucar.edu" target="_blank">brownrig@ucar.edu</a>>:<br>
<br>
> Hi Stavros,<br>
><br>
> To the best of my knowledge, the masking functions are unaware of any<br>
> cartographic projections in the data, and so both grid and shapefile need<br>
> to be in the same coordinate system.  I would guess your shapefile is in<br>
> lat/lon?  NCL has an unadvertised interface to the well-known projection<br>
> library Proj4 -- if you know the projection parameters of your<br>
> stereographic grid, you could use the proj4 interface to convert the grid<br>
> to lat/lon. The function is:<br>
><br>
>   transform_coordinate(srcProj:<wbr>string, dstProj:string, x:double, y:double,<br>
> z:double)<br>
><br>
> where the src/dstProj are proj4 descriptors of the source/destination<br>
> projections. The coordinates are passed in and transformed out through the<br>
> x,y,z arrays (I can get you started if you send me the grid's metadata<br>
> describing the projection).<br>
><br>
> As for speeding up the process, one strategy I've used in the past is, for<br>
> coastal nations that have many small islands (relative to the grid<br>
> resolution), exclude those. Shapefiles have an "area" variable, and<br>
> magnitude of that value will be vastly larger than those for islands. So by<br>
> inspecting the data, you get a feel for which shapefile polygons are the<br>
> ones "of interest", and you can select those out, and then use the<br>
> gc_inout() function to perform the masking.  Alternatively, if you have<br>
> access to GIS software like ArcGIS, you could do the same culling<br>
> operations (likely much easier!) to create such a paired-down shapefile,<br>
> and then use shapefile_mask_data() as before.<br>
><br>
> Hope that helps...<br>
> Rick<br>
><br>
> On Wed, Oct 18, 2017 at 4:33 AM, Stavros Dafis <<a href="mailto:sdafis@cc.uoi.gr" target="_blank">sdafis@cc.uoi.gr</a>> wrote:<br>
><br>
> > Dear NCL users,<br>
> ><br>
> > I have a set of netcdf files with snow, land and sea variables and the<br>
> > projection is Polar Stereographic ellipsoidal (WGS-84 ellipsoid). I have<br>
> > plotted an example remapped to a Mercator projection:<br>
> > <a href="https://ibb.co/mFZAd6" rel="noreferrer" target="_blank">https://ibb.co/mFZAd6</a><br>
> ><br>
> > I would like to mask the data only inside the shapefile of Greece, but I<br>
> > cannot<br>
> > figure out why the "shapefile_mask_data" keeps the wrong values. This is<br>
> > what I<br>
> > get:<br>
> > <a href="https://ibb.co/iyBZ5m" rel="noreferrer" target="_blank">https://ibb.co/iyBZ5m</a><br>
> ><br>
> > How is it possible to cut the values in another projection than the<br>
> > initial?<br>
> ><br>
> > This is part of the code for masking:<br>
> ><br>
> > shp_filename1 = "/data/dafis/Greek_shape/GRC_a<wbr>dm1.shp"    ; State outlines<br>
> ><br>
> >   opt             = True<br>
> >   ;opt@delta_kilometers  = 2<br>
> >   ;opt@minlat      = minlat<br>
> >   ;opt@maxlat      = maxlat<br>
> >   ;opt@minlon      = minlon<br>
> >   ;opt@maxlon      = maxlon<br>
> >   opt@debug       = True<br>
> >   opt@shape_var   = "NAME_0"<br>
> >   opt@shape_names = "Greece"<br>
> >   ;opt@return_mask = True<br>
> >   opt@keep        = True<br>
> ><br>
> >  snow_mask = shapefile_mask_data(snow,shp_f<wbr>ilename1,opt)<br>
> ><br>
> > The whole procedure takes 5 minutes for one day and I have 6205 days to<br>
> > compute.<br>
> > I read in this thread that is normal to take so long to run, due to the<br>
> > complex<br>
> > coastline:<br>
> > <a href="https://www.ncl.ucar.edu/Support/talk_archives/2014/0728.html" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Suppo<wbr>rt/talk_archives/2014/0728.<wbr>html</a><br>
> ><br>
> > Any suggestions? The opt@minlat/lon etc does not help to plot faster.<br>
> ><br>
> > Has anyone succeeded in speeding up the procedure or do you know where I<br>
> > can<br>
> > find coarse shapefiles for countries since I mostly care about excluding<br>
> > the<br>
> > neighbor countries and not the sea.<br>
> ><br>
> ><br>
> > --<br>
> > Stavros NTAFIS (DAFIS)<br>
> > ------------------------------<wbr>-----------------<br>
> > Physicist - Meteorologist, M.Sc.<br>
> > Ph.D. Candidate, Polytechnic School of Paris, France<br>
> > Laboratory of Dynamic Meteorology (LMD)<br>
> > Research Associate, National Observatory of Athens<br>
> > Tel. : <a href="tel:%28%2B33%299%2081%2094%2022%2012" value="+33981942212" target="_blank">(+33)9 81 94 22 12</a><br>
> > Mobile: <a href="tel:%28%2B33%29%20066%20030%200147" value="+33660300147" target="_blank">(+33) 066 030 0147</a><br>
> > <a href="tel:%28%2B30%29%20697%2004%2020%20242" value="+306970420242" target="_blank">(+30) 697 04 20 242</a><br>
> > ---------------<br>
> > Weather charts:<br>
> > <a href="http://www.meteo.gr" rel="noreferrer" target="_blank">http://www.meteo.gr</a><br>
> > <a href="http://www.meteo.gr/meteomaps/" rel="noreferrer" target="_blank">http://www.meteo.gr/meteomaps/</a><br>
> ><br>
> ><br>
> ><br>
> > ______________________________<wbr>_________________<br>
> > ncl-talk mailing list<br>
> > <a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
> > List instructions, subscriber options, unsubscribe:<br>
> > <a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailma<wbr>n/listinfo/ncl-talk</a><br>
> ><br>
><br>
</blockquote></div><br></div></div></div></div></div></div></div></div></div></div>