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

Rick Brownrigg brownrig at ucar.edu
Mon Nov 13 08:33:50 MST 2017 

Hi Stavros,

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!

Rick

On Sun, Nov 12, 2017 at 12:32 PM, Dennis Shea <shea at ucar.edu> wrote:

> FYI: netCDF supports the following types:
>
> http://www.unidata.ucar.edu/software/netcdf/docs/data_type.html
>
> ---
> It dies not support type "logical".
>
> I suggest a simple integer (eg: type byte if file size is important)
>
> Your current code:
>
> snow_mask = new(dimsizes(snow), logical)
> snow_mask =False
>
> do i=0,689
>   startI = b->segments(i,0)
>   endI = startI + b->segments(i,1) - 1
>   if ((endI-startI).lt.10000) then
>     continue
>   end if
>   print("area: " + i)
>   snow_mask = snow_mask .or. gc_inout(lat, lon, y(startI:endI),
> x(startI:endI))
> end do
> printVarSummary(snow)
> printMinMax(snow,True)
>
> ====
>
> SNOW_MASK= where(snow_mask, tobyte(1), tobyte(0))
> copy_VarCoords(snow_mask, SNOW_MASK)
> SNOW_MASK at long_name = snow_mask at long_name
> SNOW_MASK at info   = "1=True; 0=False"
>
>
>
>
>
>
>
>
> On Sun, Nov 12, 2017 at 3:37 AM, Stavros Dafis <sdafis at cc.uoi.gr> wrote:
>
>> 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/AACGM3BfP26G9VgMU
>> QeuD1kDa?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
>> > > > >
>> > > >
>> > >
>> >
>> _______________________________________________
>> 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/20171113/8a0af925/attachment.html>

More information about the ncl-talk mailing list