[ncl-talk] area of a shapefile
Mary Haley
haley at ucar.edu
Wed Mar 14 21:24:16 MDT 2018
Hi Vanúcia,
The extra shapefile files were helpful, thanks.
I've attached a heavily modified version of the script you provided to
Dennis. I think he had made some modifications before he handed it off to
me.
I tried to clean this code up and create several functions and procedures
so it is easier to see what the main
driver
code is doing. Search for the section that say
s:
This driver code does the following
:
which is the main code that calls the various functions and procedures. The
script should be self-explanatory.
I tried to approximate the area of your masked data by creating
quadrilaterals that surrounds each of the lat/lon points where the data was
non-missing. I then used the gc_qarea function to calculate the area of a
quadrilateral on a sphere. The area_poly_sphere function was used to
calculate the area of the whole shapefile.
You will need to decide if the quadrilateral method is accurate enough for
you. Note that some of the rectangles fall outside the area of the
shapefile outline (see the area_shapefile_mask_example.000003.png image).
You can fix this kind of thing by making the rectangle smaller, based on
whether it has missing values on any side of it. It's up to you how
accurate you want this to be.
I would like to make this an example on our shapefile page. Would this be
okay with you?
--Mary
On Wed, Mar 14, 2018 at 9:04 AM, Vanúcia Schumacher <
vanucia-schumacher at hotmail.com> wrote:
> Hi Mary,
>
>
> The print_shapefile_info:
>
>
> (0) ======================================================================
> (0) Filename: "Cuenca_Universidad.shp"
> (0) Geometry type: polygon
> (0) # of features: 1
> (0) Min/max lat: -34.77/ -34.62
> (0) Min/max lon: -70.44/ -70.25
> (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) ======================================================================
>
>
> The shapefiles are attached.
>
>
> 2) The inventory.asc it's a mascara file, in which values are only for
> some pixels of interest.
>
> I need to select the non-missing values of this file for the shapefile
> region.
>
> And then calculate the area of the file, according to the original file:
> inventory.asc
>
> I want to know the value of the area for non-missing data in the shapefile
> region
> ------------------------------
> *De:* Mary Haley <haley at ucar.edu>
> *Enviado:* quarta-feira, 14 de março de 2018 11:42:58
> *Para:* Vanúcia Schumacher
> *Cc:* Dennis Shea; ncl-talk at ucar.edu
>
> *Assunto:* Re: [ncl-talk] area of a shapefile
>
> Vanúcia,
>
> Dennis asked me to look into this.
>
> I tried looking at the shapefile and script you provided, but either
> something is wrong with the shapefile, or we don't have all the necessary
> files. Without the necessary files, you cannot correctly mask your data
> against a shapefile.
>
> I used a simple procedure called "print_shapefile_info" to print out basic
> information about the shapefile, and the lat/lon values are out-of-range. I
> think this is because the Cuenca_Universidad.prj file is not available.
> Can you provide this?
>
> Here is the simple script I used, by downloading "shapefiles_utils.ncl"
> from the shapefiles example page:
>
> load "./shapefile_utils.ncl"
> sname = "Cuenca_Universidad.shp"
> print_shapefile_info(sname)
> ; plot_shapefile(sname)
>
> And here's the output. Note the min/max of the lat/lons, which are very
> clearly wrong. They should be in the range -90 to 90 and -180 to 180:
>
> (0) Filename: "Cuenca_Universidad.shp"
> (0) Geometry type: polygon
> (0) # of features: 1
> (0) Min/max lat: 6151442.45/6167964.95
> (0) Min/max lon: 368312.46/385366.64
> (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
>
> Secondly, Dennis and I are not sure what you are trying to accomplish with
> the script. When the inventory values are read in, many of them are
> missing, and the rest appear to be all equal to 1. Can you explain more
> about what you want to do with this data?
>
> Thanks,
>
>
> --Mary
>
>
>
> On Mon, Mar 12, 2018 at 12:50 PM, Vanúcia Schumacher <
> vanucia-schumacher at hotmail.com> wrote:
>
> Thanks for the tips, but I'm not making progress on this issue.
>
>
> I tried to follow example 4 of the function: area_poly_sphere
> <https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml>
> ,
>
> but the area I need to calculate is not from the shapefile, it refers to inventory.asc
> file,
> in which I use the shapefile to extract the pixel values in the polygon
> region of the shapefile.
> The only information I have from the inventory.asc are:
>
> ncols 129
> nrows 156
> xllcorner -70.8944702148438
> yllcorner -35.3186645507812
> xcellsize 0.0109375946281492
> ycellsize 0.00894945095746945
> nodata_value -9999
>
> The script I tried follows below, but I get the following error
> gc_clkwise: the arrays must have the same dimension sizes
>
> link to line:
> if (gc_clkwise(mrb_lat, mrb_lon)) then
>
> which I believe that lat and lon information are wrong.
>
> Thank you very much any help
>
>
>
> Variable: lonss
> Type: float
> Total Size: 516 bytes
> 129 values
> Number of Dimensions: 1
> Dimensions and sizes: [129]
> Coordinates:
>
> Variable: latss
> Type: float
> Total Size: 624 bytes
> 156 values
> Number of Dimensions: 1
> Dimensions and sizes: [156]
> Coordinates:
> (0) n_mrb=156
>
>
> Script:
>
> ncols = 129
> nrows = 156
> lonLL = -70.8944702148438
> latLL = -35.3186645507812
> deltaLon = 0.0109375946281492
> deltaLat = 0.00894945095746945
> missingVal = -9999.0
>
> data = readAsciiTable("inventory.asc",ncols,"float",7)
>
> ;----- construct 1D coordinates....
> lons = ispan(0,ncols-1,1) * deltaLon + lonLL
> lats = ispan(0,nrows-1,1) * deltaLat + latLL
> lats = lats(::-1)
>
> ;---- gsn_coordinates wants 2D coordinate arrays
> data at lon2d = conform_dims((/ nrows, ncols/), lons, 1)
> data at lat2d = conform_dims((/ nrows, ncols/), lats, 0)
> data at _FillValue = missingVal
>
> ;---- Open shapefile .
>
> shp_filename = "universidad.shp"
> data_mask = shapefile_mask_data(data,shp_filename,True)
>
> ;---- Selecting just non-missing values, from the inventory.asc using
> shapefile region:
>
> x = ndtooned(data_mask(:,:)) ; reduce to 1D array
> igood = ind(.not.ismissing(x)) ; extract only the
> non-missing values.
> xgood = x(igood)
>
> ;----- construct 1D coordinates to xgood var:
>
> lonss = ispan(0,ncols-1,1) * deltaLon + lonLL
> latss = ispan(0,nrows-1,1) * deltaLat + latLL
> latss = lats(::-1)
>
> ;---- gsn_coordinates wants 2D coordinate arrays to xgood var:
>
> xgood at lon2d = conform_dims((/ nrows, ncols/), lons, 1)
> xgood at lat2d = conform_dims((/ nrows, ncols/), lats, 0)
> xgood at _FillValue = missingVal
>
>
> mrb_lon = lonss Problem is here!
> mrb_lat = latss
>
> n_mrb = dimsizes(mrb_lat)
> print("n_mrb="+n_mrb)
>
> re = 6371.
> re at units= "km"
>
> if (gc_clkwise(mrb_lat, mrb_lon)) then
> area_mrb = area_poly_sphere(mrb_lat, mrb_lon, re)
>
> area_mrb at units = "km^2"
> area_mrb at long_name = "Area Enclosed by the MRB"
>
> print (area_mrb)
> else
> print ("MRB points are not in clockwise order.")
> exit
> end if
>
>
>
>
> <https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml>
>
> ------------------------------
> *De:* Dennis Shea <shea at ucar.edu>
> *Enviado:* segunda-feira, 12 de março de 2018 11:25:03
> *Para:* Vanúcia Schumacher
> *Cc:* ncl-talk at ucar.edu
> *Assunto:* Re: [ncl-talk] area of a shapefile
>
> You could click "Functions" [Alphabetical]; then, search for the word
> 'area'
> Keep clicking ....
>
> https://www.ncl.ucar.edu/Document/Functions/Built-in/area_
> poly_sphere.shtml
>
> Also,
> https://www.ncl.ucar.edu/Document/Functions/Built-in/gc_qarea.shtml
> See Examples
>
> Good luck
>
>
> On Mon, Mar 12, 2018 at 8:14 AM, Vanúcia Schumacher <
> vanucia-schumacher at hotmail.com> wrote:
>
> Hi NCL users,
>
> I would like help to calculate the area from the shapefile.
> Any ideas how can I proceed?
>
>
> hdr = readAsciiHead("inventory.asc", 7)
> ncols = stringtoint( str_get_field(hdr(0), 2, " ") )
> nrows = stringtoint( str_get_field(hdr(1), 2, " ") )
> lonLL = stringtofloat( str_get_field(hdr(2), 2, " ") )
> latLL = stringtofloat( str_get_field(hdr(3), 2, " ") )
> deltaLon = stringtofloat( str_get_field(hdr(4), 2, " " ) )
> deltaLat = stringtofloat( str_get_field(hdr(5), 2, " ") )
> missingVal = stringtofloat( str_get_field(hdr(6), 2, " ") )
>
> data = readAsciiTable("inventory.asc",ncols,"float",7)
> data at _FillValue =-9999.0
>
> ;----- construct 1D coordinates....
> lons = ispan(0,ncols-1,1) * deltaLon + lonLL
> lats = ispan(0,nrows-1,1) * deltaLat + latLL
> lats = lats(::-1)
>
> ;---- gsn_coordinates wants 2D coordinate arrays
> data at lon2d = conform_dims((/ nrows, ncols/), lons, 1)
> data at lat2d = conform_dims((/ nrows, ncols/), lats, 0)
> data at _FillValue = missingVal
>
> ;---- Open shapefile and read lat/lon values.
>
> shp_filename = "Cuenca_Universidad.shp"
> data_mask = shapefile_mask_data(data,shp_filename,True)
>
> ;--- Area ?
>
>
> _______________________________________________
> 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/20180314/b4a8823a/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: area_shapefile_mask_example.000001.png
Type: image/png
Size: 41333 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/b4a8823a/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: area_shapefile_mask_example.000002.png
Type: image/png
Size: 72788 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/b4a8823a/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: area_shapefile_mask_example.000003.png
Type: image/png
Size: 92079 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/b4a8823a/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: area.shapefile2_mod.ncl
Type: application/octet-stream
Size: 12302 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/b4a8823a/attachment.obj>
More information about the ncl-talk
mailing list