[ncl-talk] area of a shapefile
Vanúcia Schumacher
vanucia-schumacher at hotmail.com
Thu Mar 15 06:03:24 MDT 2018
Thank you so much Mary and Dennis for the support.
Of course, you can make this example on the page.
________________________________
De: Mary Haley <haley at ucar.edu>
Enviado: quinta-feira, 15 de março de 2018 00:24:16
Para: Vanúcia Schumacher
Cc: ncl-talk at ucar.edu
Assunto: Re: [ncl-talk] area of a shapefile
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<mailto: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<mailto: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<mailto: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/20180315/2a47a199/attachment-0001.html>
More information about the ncl-talk
mailing list