[ncl-talk] Masking within a shapefile Polar Stereographic projection
Dennis Shea
shea at ucar.edu
Sun Nov 12 12:32:13 MST 2017
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/
> 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
> > > > >
> > > >
> > >
> >
> _______________________________________________
> 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/20171112/343995d4/attachment.html>
More information about the ncl-talk
mailing list