[ncl-talk] shapefiles

Dave Allured - NOAA Affiliate dave.allured at noaa.gov
Fri Feb 8 20:13:09 MST 2019


Vanúcia,

Example shapefiles_21 does not work the way I expected.  It needs a
shapefile with areas surrounding the region of interest, because it works
by blanking out unwanted areas with white, leaving only the central area
visible.  NCL does not currently have the ability to show contour plot
areas inside specified polygons.  It can only blank out unwanted polygons
in graphics mode.

The shapefile in website example shapefiles_21 contains surrounding areas,
as seen in figure 1 below.  When you specify the list of area names to
display, what actually happens is that all of the areas NOT in the list are
blanked out in white.  Please notice the difference between this shapefile,
and your original shapefile that contains only your wanted areas.

Your original program crashed because in "areas_of_interest" you listed ALL
AREAS inside the shape file.  This reduced the number of blanking areas to
zero.  This resulted in an unexpected zero length array which was not
properly handled by the program.  There is nothing basically wrong with the
shapefile, it's just that example shapefiles_21 can't use it in this manner.

To demonstrate, I removed one of the four areas from your area list.  This
avoided the program crash.  The modified script with a few other changes is
attached.  In figure 2, notice two things.  (1) The area Cuenca Coton is
blanked to white, because I removed it from the area list.  (2) All of the
unwanted surrounding areas are VISIBLE in the contour plot, because they
are NOT available in the shapefile; therefore they can not be blanked out.
I hope this shows you how example 21 actually works.

(A) Your best solution would be to find a larger shapefile which contains
many basin outlines for this part of Chile.  That should just drop in to
the current script, and blank out the surrounding unwanted areas.

(B) It is possible to construct exclusion masks using shapefile editing
tools.  Draw two rectangles which completely cover your plot frame, then
use your original shapefile to cut out the desired visible areas, leaving
only exclusion areas.  You need two rectangles, not one, because NCL does
not currently handle polygons with internal holes.  There are many
shapefile tools, but I don't have experience with them.  Perhaps someone
else could recommend a tool that can do this.

(C) It is even possible to use NCL to manually extract polyline segments
from your original shapefile, and combine them with rectangle vertices to
construct these two exclusion masks.  This is tricky and tedious, so the
first two solutions are much better if you can do them.  Good luck.

--Dave


On Thu, Jan 31, 2019 at 5:58 AM Vanúcia Schumacher <
vanucia-schumacher at hotmail.com> wrote:

> Unfortunately I do not have another shapefile, just this one referring to
> my work.
>
> *De:* Dave Allured - NOAA Affiliate <dave.allured at noaa.gov>
> *Enviado:* terça-feira, 29 de janeiro de 2019 16:16
>
> Vanúcia,
>
> If you can find another shape file of the same region from a different
> source, give that a try.  This might be quicker than waiting for analysis
> of the original shape file.  If you find a working shape file, please let
> us know.
>
> --Dave
>
>
> On Mon, Jan 28, 2019 at 10:50 AM Dave Allured - NOAA Affiliate <
> dave.allured at noaa.gov> wrote:
>
> Thanks.  I can't look at this today.  I will get back to you in a day or
> two.  Perhaps someone else would like to check out the shape file and the
> script.
>
>
> On Mon, Jan 28, 2019 at 10:45 AM Vanúcia Schumacher <
> vanucia-schumacher at hotmail.com> wrote:
>
> Is attached.
>
> Thank you for the advance
>
> *De:* Dave Allured - NOAA Affiliate <dave.allured at noaa.gov>
> *Enviado:* segunda-feira, 28 de janeiro de 2019 15:37
>
> Vanúcia,
>
> Shape files come in different varieties.  It looks like example
> shapefiles_21 may need an upgrade to be able to read your shape file.  If
> you can send a copy of the shape file, and also your current script
> shape21.ncl, I will test it on my system.  Please do not include the
> scientific data file; I do not need it for testing the shape file.  The NCL
> support team might need to fix this.
>
> --Dave
>
>
> On Mon, Jan 28, 2019 at 9:47 AM Vanúcia Schumacher <
> vanucia-schumacher at hotmail.com> wrote:
>
> Thanks for the tip.
> I'm trying to use the script from example shapefiles_21, but I'm having
> this error in function. I can not fix it.
>
> fatal:Subscript out of range, error in subscript #0
> fatal:An error occurred reading features
> fatal:["Execute.c":8640]:Execute: Error occurred at or near line 56 in
> file shape21.ncl
> fatal:["Execute.c":8640]:Execute: Error occurred at or near line 250 in
> file shape21.ncl
>
> The lines correspond to:
>
> function get_areas_of_interest(shp_file_name,shp_var_name,opt[1]:logical)
> begin
> ;---Open the shapefile
>   f = addfile(shp_file_name,"r")
>   features = f->$shp_var_name$
>
>   if(opt.and.isatt(opt,"areas_to_exclude")) then
>     features at _FillValue = default_fillvalue(typeof(features))
>     do na=0,dimsizes(opt at areas_to_exclude)-1
>       ii := ind(features.eq.opt at areas_to_exclude(na))
>       if(.not.any(ismissing(ii))) then
>         features(ii) = features at _FillValue
>       end if
>     end do
>     return(features(ind(.not.ismissing(features))))    ------>  line 56
>   else
>     return(features)
>   end if
> end
>
> opt                  = True
>  opt at areas_to_exclude = areas_of_interest
>   areas_to_fill        =
> get_areas_of_interest(shp_filename1,shape_var_name,opt)    ----> line 250
>
> Infos:
>
> 0) Terrain Height (m) : min=0   max=5070.39
> (0) ======================================================================
> (0) Filename: "cuencasx.shp"
> (0)    Geometry type: polygon
> (0)    # of features: 4
> (0)    Min/max lat:    -34.77/ -34.23
> (0)    Min/max lon:    -70.53/ -70.13
> (0)    Variable names and their types:
> (0)        geometry : integer
> (0)        segments : integer
> (0)        x : double
> (0)        y : double
> (0)        COD_CUEN : string
> (0)        COD_SUBC : string
> (0)        COD_SSUBC : string
> (0)        NOMBRE : string
> (0)        AREAKM2 : double
> (0) ======================================================================
>
> Variable: features
> Type: string
> Total Size: 32 bytes
>             4 values
> Number of Dimensions: 1
> Dimensions and sizes: [num_features | 4]
> Coordinates:
> Number Of Attributes: 0
> (0) Cuenca Cipreses
> (1) Cuenca Coton
> (2) Cuenca Cortaderal
> (3) Cuenca Universidad
> (0) ==================================================
> (0) Shapefile:         cuencasx.shp
> (0) Areas of interest: Cuenca Cipreses,Cuenca Coton,Cuenca
> Cortaderal,Cuenca Universidad
> (0) min_lat_chk:       -37.4943
> (0) max_lat_chk:       -31.9474
> (0) min_lon_chk:       -73.6526
> (0) max_lon_chk:       -66.7254
> (0) min_lat_data:      -37.4943
> (0) max_lat_data:      -31.9474
> (0) min_lon_data:      -73.6526
> (0) max_lon_data:      -66.7254
> (0) 4761 data values originally
> (0) Will keep data values inside given shapefile areas
> (0) ==================================================
> (0) 13 data values kept
> (0) shapefile_mask_data: elapsed time: 0.00758901 CPU seconds.
> (0) ==================================================
>
>
> *De:* Dave Allured - NOAA Affiliate <dave.allured at noaa.gov>
> *Enviado:* domingo, 27 de janeiro de 2019 18:23
>
> Vanúcia,
>
> Your two data files are on different grids.  If you only need to make
> plots that show the same region on a map, then use graphics masking, not
> data masking.  See the third plot in example shapefiles_21 on this page:
>
> https://www.ncl.ucar.edu/Applications/shapefiles.shtml
>
> If you need to perform calculations between the two data sets, then you
> will need to regrid the data from one file to match the other grid.  Please
> see this documentation for regridding:
>
> https://www.ncl.ucar.edu/Applications/regrid.shtml
> https://www.ncl.ucar.edu/Applications/ESMF.shtml
>
> --Dave
>
>
> On Sun, Jan 27, 2019 at 12:40 PM Vanúcia Schumacher <
> vanucia-schumacher at hotmail.com> wrote:
>
> Hi all,
>
> I need support to find the problem with my script, with the purpose of
> cutting two different data in a shapefile, but keeping the information
> such as time lat and lon.
> My script (attachment) is "cutting" the region from shapefile different
> for each input (data) (see Figure attachment).
> I'd like to select the shapefile region of the same size independent of
> the input file, and keep the time, lat, and lon information.
>
> I appreciate any help
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190208/d7991025/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gadm36_USA_1.png
Type: image/png
Size: 88175 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190208/d7991025/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: shape21c.000003.png
Type: image/png
Size: 178440 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190208/d7991025/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: shape21c.ncl
Type: application/octet-stream
Size: 11767 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190208/d7991025/attachment-0001.obj>


More information about the ncl-talk mailing list