[ncl-talk] area of a shapefile

[Vanúcia Schumacher]

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<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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180312/ba3e70ab/attachment.html>