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

Wed Oct 18 12:41:04 MDT 2017

Stavros and Rick,

First, to address the slowness: in addition to what Rick suggested: if the
masking locations don't change for the 6000+ time steps, then what you can
do is have the shapefile mask function return a mask array of 0s and 1s,
and then you can simply apply this mask using the "where" function, instead
of calling shapefile_mask_data every time. This is what shapefile_20.ncl is
doing:

http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex20

With regard to the image being wrong, it looks more like a plotting issue
than a masking issue. It looks like either the lat/lon arrays are not being
passed to the plotting function, OR, you are trying to do a "native"
projection by setting tfDoNDCOverlay to True, and then setting the map
projection parameters yourself, but the projection is not correct.

My guess is that "snow_mask" doesn't have the necessary lat/lon data
associated with it, and you may need to copy it over from "snow".

Can you provide a script so we can see what kind of plotting you're doing?

Thanks,

--Mary

On Wed, Oct 18, 2017 at 7:49 AM, Rick Brownrigg <brownrig at ucar.edu> wrote:

> 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
>>
>
>
> _______________________________________________
> 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/9025e7b3/attachment.html>