[ncl-talk] area of a shapefile
Mary Haley
haley at ucar.edu
Wed Mar 14 08:42:58 MDT 2018
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/7c0784f9/attachment.html>
More information about the ncl-talk
mailing list