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

Rick Brownrigg brownrig at ucar.edu
Fri Nov 3 11:11:24 MDT 2017 

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
> > >
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171103/46f8131e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Mask_snow.ncl
Type: application/octet-stream
Size: 2026 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171103/46f8131e/attachment.obj>

More information about the ncl-talk mailing list