[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