<div dir="ltr"><div><div>Hi Stavros,<br><br></div>My apologies -- I had not tested my suggestion end-to-end, and I was unaware of the NetCDF limitation. Your solution looks like a good one however, and thanks for contributing that back to the group!<br><br></div>Rick<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Nov 12, 2017 at 12:32 PM, Dennis Shea <span dir="ltr"><<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div>FYI: netCDF supports the following types:<br><br><a href="http://www.unidata.ucar.edu/software/netcdf/docs/data_type.html" target="_blank">http://www.unidata.ucar.edu/<wbr>software/netcdf/docs/data_<wbr>type.html</a><br><br>---<br></div>It dies not support type "logical".<br><br></div>I suggest a simple integer (eg: type byte if file size is important)<br><br>Your current code:<br><br>snow_mask = new(dimsizes(snow), logical)<br>snow_mask =False<br><br>do i=0,689<br> startI = b->segments(i,0)<br> endI = startI + b->segments(i,1) - 1<br> if ((endI-startI).lt.10000) then<br> continue<br> end if<br> print("area: " + i)<br> snow_mask = snow_mask .or. gc_inout(lat, lon, y(startI:endI), x(startI:endI))<br>end do<br>printVarSummary(snow)<br>printMinMax(snow,True)<br><br>====<br><br>SNOW_MASK= where(snow_mask, tobyte(1), tobyte(0))<br>copy_VarCoords(snow_mask, SNOW_MASK)<br></div>SNOW_MASK@long_name = snow_mask@long_name<br></div>SNOW_MASK@info = "1=True; 0=False"<br><br><br><div><div><br><br><br><div><div><br><br></div></div></div></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Nov 12, 2017 at 3:37 AM, 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:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Dear Rick,<br>
<br>
thank you again for your help. Just one last question about writing in netCDF. I<br>
often use to write netCDF files with NCL but this is my first time trying to<br>
write a logical (character) array (False/True). I found a thread back in 2014<br>
about this problem:<br>
<a href="https://www.ncl.ucar.edu/Support/talk_archives/2014/1003.html" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Suppo<wbr>rt/talk_archives/2014/1003.<wbr>html</a><br>
<br>
Is it solved in NCL 6.2.0 which I am using or is it still a bug?<br>
<br>
regards,<br>
<span class="m_-7537011226744105384HOEnZb"><font color="#888888"><br>
--<br>
Stavros<br>
</font></span><div class="m_-7537011226744105384HOEnZb"><div class="m_-7537011226744105384h5"><br>
<br>
<br>
Quoting Rick Brownrigg <<a href="mailto:brownrig@ucar.edu" target="_blank">brownrig@ucar.edu</a>>:<br>
<br>
> Hi Stavros,<br>
><br>
> I'm replying back to the list as the technique I'm suggesting may be useful<br>
> to others.<br>
><br>
> As I understand your original problem was that masking was taking a long<br>
> time, and you need to apply that make to a large number of datasets. I<br>
> found that with the shapefile you sent, it was taking hours, not minutes. A<br>
> simpler, more generalized shapefile would go a long ways towards<br>
> alleviating the problem. But, given the shapefile at hand...<br>
><br>
> My thought was to create the mask once, and write it to a NetCDF file. Then<br>
> when processing the individual files, load that mask and apply it against<br>
> the data. The attached script shows most of the steps involved (it does not<br>
> write to NetCDF, but that's simple enough -- examples at<br>
> <a href="http://ncl.ucar.edu/Applications/o-netcdf.shtml" rel="noreferrer" target="_blank">http://ncl.ucar.edu/Applicatio<wbr>ns/o-netcdf.shtml</a>).<br>
><br>
> - reads the lat/lon and snow data. Gets the x/y from the shapefile.<br>
> - creates a boolean "mask" variable the same dimensions as the snow data<br>
> - Loops over each of the individual polygons from the shapefile, and calls<br>
> gc_inout() to find what part of the data grid overlaps individual polygons.<br>
> gc_inout() returns a logical array denoting which grid cells are in/out.<br>
> Those results are logically OR'd into the mask. On the premise that the<br>
> smaller islands are not of concern relative to the grid resolution, the<br>
> in/out test is skipped for polygons that have fewer vertices than a<br>
> threshold (arbitrarily, 10,000 here -- that dragged in only the mainland<br>
> and 3 other larger islands). Even so, this loop/mask-generation phase<br>
> takes 2-3hrs with the given shapefile.<br>
> - 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>
> This assumes your data don't change spatial extent from file to file. If<br>
> this approach is useful, you'd want to break the script into two pieces: i)<br>
> to generate and save the mask, ii) read your data, apply the mask, and<br>
> generate the plots.<br>
><br>
> BTW: as a side comment, I may have thrown noise into the previous<br>
> discussion regarding map projections. Your data and shapefile were both in<br>
> geographic coordinates.<br>
><br>
> Hope this helps...<br>
> Rick<br>
><br>
> On Wed, Nov 1, 2017 at 12:27 PM, Stavros Dafis <<a href="mailto:sdafis@cc.uoi.gr" target="_blank">sdafis@cc.uoi.gr</a>> wrote:<br>
><br>
> > Hello Rick,<br>
> ><br>
> > Thank you for answering my question about masking data in NCL. I am sorry<br>
> > 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<br>
> > 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,<br>
> > 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<br>
> > 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,<br>
> > 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<br>
> > 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<br>
> > 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<br>
> > I<br>
> > > > cannot<br>
> > > > figure out why the "shapefile_mask_data" keeps the wrong values. This<br>
> > 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<br>
> > 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<br>
> > I<br>
> > > > can<br>
> > > > find coarse shapefiles for countries since I mostly care about<br>
> > 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>
> ><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>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>