[ncl-talk] area of a shapefile
Vanúcia Schumacher
vanucia-schumacher at hotmail.com
Wed Mar 14 09:04:02 MDT 2018
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<mailto: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<mailto: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<mailto: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<mailto: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<mailto: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<mailto: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/439b5422/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.cpg
Type: application/octet-stream
Size: 5 bytes
Desc: Cuenca_Universidad.cpg
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.dbf
Type: application/x-dbf
Size: 346 bytes
Desc: Cuenca_Universidad.dbf
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.prj
Type: application/octet-stream
Size: 409 bytes
Desc: Cuenca_Universidad.prj
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.sbn
Type: application/octet-stream
Size: 132 bytes
Desc: Cuenca_Universidad.sbn
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.sbx
Type: application/octet-stream
Size: 116 bytes
Desc: Cuenca_Universidad.sbx
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment-0003.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.shp
Type: application/octet-stream
Size: 3372 bytes
Desc: Cuenca_Universidad.shp
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment-0004.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Cuenca_Universidad.shx
Type: application/octet-stream
Size: 108 bytes
Desc: Cuenca_Universidad.shx
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/439b5422/attachment-0005.obj>
More information about the ncl-talk
mailing list