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

[Rick Brownrigg]

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/20171018/8496602c/attachment.html>