[ncl-talk] Masking within a shapefile Polar Stereographic projection

Stavros Dafis sdafis at cc.uoi.gr
Sun Nov 12 03:37:17 MST 2017 

Dear Rick,

thank you again for your help. Just one last question about writing in netCDF. I
often use to write netCDF files with NCL but this is my first time trying to
write a logical (character) array (False/True). I found a thread back in 2014
about this problem:
https://www.ncl.ucar.edu/Support/talk_archives/2014/1003.html

Is it solved in NCL 6.2.0 which I am using or is it still a bug?

regards,

-- 
Stavros



Quoting Rick Brownrigg <brownrig at ucar.edu>:

> Hi Stavros,
>
> I'm replying back to the list as the technique I'm suggesting may be useful
> to others.
>
> 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...
>
> 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
> http://ncl.ucar.edu/Applications/o-netcdf.shtml).
>
> - reads the lat/lon and snow data. Gets the x/y from the shapefile.
> - creates a boolean "mask" variable the same dimensions as the snow data
> - 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.
> - Use the mask and the where() function to mask out the data:
>
>     snow = where(snow_mask.eq.True, snow, snow at _FillValue)
>
> 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.
>
> 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.
>
> Hope this helps...
> Rick
>
> On Wed, Nov 1, 2017 at 12:27 PM, Stavros Dafis <sdafis at cc.uoi.gr> wrote:
>
> > Hello Rick,
> >
> > Thank you for answering my question about masking data in NCL. I am sorry
> > for
> > not replying earlier but I experienced some problems with my machines and I
> > wanted to install Proj.4 and try your solution.
> >
> > I have not succeeded in converting the projection of my data to the
> > projection
> > of the shapefile. Maybe I am using Proj.4 in a wrong way and I found no
> > particular documentation about the function "transform_coordinate".
> >
> > My data have the following projection:
> >
> > "Polar stereographic ellipsoidal projection with WGS-84 ellipsoid"
> >
> > I have uploaded an example of the data in netCDF, my NCL script and the
> > shapefile of Greece:
> > https://www.dropbox.com/sh/rfp5yxk1abfmmgf/AACGM3BfP26G9VgMUQeuD1kDa?dl=0
> >
> > Any help would be very much appreciated.
> >
> > --
> > Stavros NTAFIS (DAFIS)
> > -----------------------------------------------
> > Physicist - Meteorologist, M.Sc.
> > Ph.D. Candidate, Polytechnic School of Paris, France
> > Laboratory of Dynamic Meteorology (LMD)
> > Research Associate, National Observatory of Athens (NOA/IERSD)
> > Tel. : (+33)1 69 33 51 64
> > (+33)9 81 94 22 12
> > Mobile: (+33) 066 030 0147
> > (+30) 697 04 20 242
> > ---------------
> > Weather charts:
> > http://www.meteo.gr
> > http://www.meteo.gr/meteomaps/
> >
> >
> > Quoting Rick Brownrigg <brownrig at ucar.edu>:
> >
> > > Hi Stavros,
> > >
> > > To the best of my knowledge, the masking functions are unaware of any
> > > cartographic projections in the data, and so both grid and shapefile need
> > > to be in the same coordinate system.  I would guess your shapefile is in
> > > lat/lon?  NCL has an unadvertised interface to the well-known projection
> > > library Proj4 -- if you know the projection parameters of your
> > > stereographic grid, you could use the proj4 interface to convert the grid
> > > to lat/lon. The function is:
> > >
> > >   transform_coordinate(srcProj:string, dstProj:string, x:double,
> > y:double,
> > > z:double)
> > >
> > > where the src/dstProj are proj4 descriptors of the source/destination
> > > projections. The coordinates are passed in and transformed out through
> > the
> > > x,y,z arrays (I can get you started if you send me the grid's metadata
> > > describing the projection).
> > >
> > > As for speeding up the process, one strategy I've used in the past is,
> > for
> > > coastal nations that have many small islands (relative to the grid
> > > resolution), exclude those. Shapefiles have an "area" variable, and
> > > magnitude of that value will be vastly larger than those for islands. So
> > by
> > > inspecting the data, you get a feel for which shapefile polygons are the
> > > ones "of interest", and you can select those out, and then use the
> > > gc_inout() function to perform the masking.  Alternatively, if you have
> > > access to GIS software like ArcGIS, you could do the same culling
> > > operations (likely much easier!) to create such a paired-down shapefile,
> > > and then use shapefile_mask_data() as before.
> > >
> > > Hope that helps...
> > > Rick
> > >
> > > On Wed, Oct 18, 2017 at 4:33 AM, Stavros Dafis <sdafis at cc.uoi.gr> wrote:
> > >
> > > > Dear NCL users,
> > > >
> > > > I have a set of netcdf files with snow, land and sea variables and the
> > > > projection is Polar Stereographic ellipsoidal (WGS-84 ellipsoid). I
> > have
> > > > plotted an example remapped to a Mercator projection:
> > > > https://ibb.co/mFZAd6
> > > >
> > > > I would like to mask the data only inside the shapefile of Greece, but
> > I
> > > > cannot
> > > > figure out why the "shapefile_mask_data" keeps the wrong values. This
> > is
> > > > what I
> > > > get:
> > > > https://ibb.co/iyBZ5m
> > > >
> > > > How is it possible to cut the values in another projection than the
> > > > initial?
> > > >
> > > > This is part of the code for masking:
> > > >
> > > > shp_filename1 = "/data/dafis/Greek_shape/GRC_adm1.shp"    ; State
> > outlines
> > > >
> > > >   opt             = True
> > > >   ;opt at delta_kilometers  = 2
> > > >   ;opt at minlat      = minlat
> > > >   ;opt at maxlat      = maxlat
> > > >   ;opt at minlon      = minlon
> > > >   ;opt at maxlon      = maxlon
> > > >   opt at debug       = True
> > > >   opt at shape_var   = "NAME_0"
> > > >   opt at shape_names = "Greece"
> > > >   ;opt at return_mask = True
> > > >   opt at keep        = True
> > > >
> > > >  snow_mask = shapefile_mask_data(snow,shp_filename1,opt)
> > > >
> > > > The whole procedure takes 5 minutes for one day and I have 6205 days to
> > > > compute.
> > > > I read in this thread that is normal to take so long to run, due to the
> > > > complex
> > > > coastline:
> > > > https://www.ncl.ucar.edu/Support/talk_archives/2014/0728.html
> > > >
> > > > Any suggestions? The opt at minlat/lon etc does not help to plot faster.
> > > >
> > > > Has anyone succeeded in speeding up the procedure or do you know where
> > I
> > > > can
> > > > find coarse shapefiles for countries since I mostly care about
> > excluding
> > > > the
> > > > neighbor countries and not the sea.
> > > >
> > > >
> > > > --
> > > > Stavros NTAFIS (DAFIS)
> > > > -----------------------------------------------
> > > > Physicist - Meteorologist, M.Sc.
> > > > Ph.D. Candidate, Polytechnic School of Paris, France
> > > > Laboratory of Dynamic Meteorology (LMD)
> > > > Research Associate, National Observatory of Athens
> > > > Tel. : (+33)9 81 94 22 12
> > > > Mobile: (+33) 066 030 0147
> > > > (+30) 697 04 20 242
> > > > ---------------
> > > > Weather charts:
> > > > http://www.meteo.gr
> > > > http://www.meteo.gr/meteomaps/
> > > >
> > > >
> > > >
> > > > _______________________________________________
> > > > ncl-talk mailing list
> > > > ncl-talk at ucar.edu
> > > > List instructions, subscriber options, unsubscribe:
> > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > > >
> > >
> >
>

More information about the ncl-talk mailing list