Re: Fwd: question about country outline

From: <brownrig_at_nyahnyahspammersnyahnyah>
Date: Mon Dec 13 2010 - 10:52:02 MST

Hi Shinn,

I took a look at the shapefile, and as we see from the
printVarSummary(f):

...
   num_features = 1 // unlimited
   num_segments = 1402
   num_points = 566678
...

This says that there's 1 *logical* polygon in the shapefile,
consisting of 1402 segments. In actuality, there are 1402 distinct
geometric polygons that comprise the Japanese islands. Unfortunately,
gc_inout() wants discrete, non-overlapping polygons, not the
shapefile's logical view. You'd in principle need to perform
gc_inout() against each of the 1402 polygons.

Depending upon the resolution of your data, you might be able to get
what you want by testing against only the largest islands. Doing the
following exploratory test:

   foo=f->segments(:,1)
   qsort(foo)
   print(foo)
...
(1394) 5793
(1395) 6603
(1396) 7663
(1397) 16566
(1398) 40122
(1399) 57991
(1400) 67561
(1401) 185495

There's a big gap between polygons(islands) with a several thousand
vertices versus those with 10,000's of vertices. Assuming those last 5
are the principle islands, you might do something like this in your
script:

   ; get islands with more than 10,000 vertices...
   mainIslands=ind(f->segments(:,1).gt.10000)

   do i=0,dimsizes(mainIslands)-1
       ; get the lon/lat coords that make up the ith geometric
polygon...
       begInd = f->segments(mainIslands(i), 0)
       endInd = begInd + f->segments(mainIslands(i), 1) - 1
       tmpX = f->x(begInd:endInd)
       tmpY = f->y(begInd:endInd)
       ...
       ... = gc_inout(...., tmpX, tmpY)
       ...
       delete(tmpX) ; in prep for next iteration
       delete(tmpY)
   end do

I hope that helps!
Rick

On Sat, 11 Dec 2010 21:35:11 +0800
  shinn wong <shinnshinnwong@gmail.com> wrote:
> ---------- Forwarded message ----------
>From: shinn wong <shinnshinnwong@gmail.com>
> Date: 2010/12/11
> Subject: Re: question about country outline
> To: Rick Brownrigg <brownrig@ucar.edu>
>
>
> Dear Rick,
> Thanks for your explanation. I've attached the shapefile (from
> http://www.diva-gis.org/gdata and then choose japan for download)
>and hope
> that you can kindly take a look at it. Below is the result of
> printVarSummary of the shapefile:
>
> Variable: f
> (0)
> filename: JPN_adm0
> path: /home/data/JPN_adm0.shp
> file global attributes:
> layer_name : JPN_adm0
> geometry_type : polygon
> geom_segIndex : 0
> geom_numSegs : 1
> segs_xyzIndex : 0
> segs_numPnts : 1
> dimensions:
> geometry = 2
> segments = 2
> num_features = 1 // unlimited
> num_segments = 1402
> num_points = 566678
> variables:
> integer geometry ( num_features, geometry )
> integer segments ( num_segments, segments )
> double x ( num_points )
> double y ( num_points )
> integer GADMID ( num_features )
> string ISO ( num_features )
> string NAME_ENGLI ( num_features )
> string NAME_ISO ( num_features )
> string NAME_FAO ( num_features )
> string NAME_LOCAL ( num_features )
> string NAME_OBSOL ( num_features )
> string NAME_VARIA ( num_features )
> string NAME_NONLA ( num_features )
> string NAME_FRENC ( num_features )
> string NAME_SPANI ( num_features )
> string NAME_RUSSI ( num_features )
> string NAME_ARABI ( num_features )
> string NAME_CHINE ( num_features )
> string WASPARTOF ( num_features )
> string CONTAINS ( num_features )
> string SOVEREIGN ( num_features )
> string ISO2 ( num_features )
> string WWW ( num_features )
> string FIPS ( num_features )
> double ISON ( num_features )
> string VALIDFR ( num_features )
> string VALIDTO ( num_features )
> double AndyID ( num_features )
> double POP2000 ( num_features )
> double SQKM ( num_features )
> double POPSQKM ( num_features )
> string UNREGION1 ( num_features )
> string UNREGION2 ( num_features )
> double DEVELOPING ( num_features )
> double CIS ( num_features )
> double Transition ( num_features )
> double OECD ( num_features )
> string WBREGION ( num_features )
> string WBINCOME ( num_features )
> string WBDEBT ( num_features )
> string WBOTHER ( num_features )
> double CEEAC ( num_features )
> double CEMAC ( num_features )
> double CEPLG ( num_features )
> double COMESA ( num_features )
> double EAC ( num_features )
> double ECOWAS ( num_features )
> double IGAD ( num_features )
> double IOC ( num_features )
> double MRU ( num_features )
> double SACU ( num_features )
> double UEMOA ( num_features )
> double UMA ( num_features )
> double PALOP ( num_features )
> double PARTA ( num_features )
> double CACM ( num_features )
> double EurAsEC ( num_features )
> double Agadir ( num_features )
> double SAARC ( num_features )
> double ASEAN ( num_features )
> double NAFTA ( num_features )
> double GCC ( num_features )
> double CSN ( num_features )
> double CARICOM ( num_features )
> double EU ( num_features )
> double CAN ( num_features )
> double ACP ( num_features )
> double Landlocked ( num_features )
> double AOSIS ( num_features )
> double SIDS ( num_features )
> double Islands ( num_features )
> double LDC ( num_features )
> double Shape_Leng ( num_features )
> double Shape_Area ( num_features )
> I've also tried to test some points. For example,
> print(gc_inout(26.7,138.4,y,x)) returns True even the point lies
>outside the
> national boundaries. Thanks once again and look forward to your
>reply.
>
> Shinn
> 2010/12/11 Rick Brownrigg <brownrig@ucar.edu>
>
> Shinn,
>>
>> Its a bit hard to say why. If I recall your previous plot
>>correctly, the
>> shapefile contained boundaries of internal regions, in addition to
>>the
>> country boundary. This likely confuses the in/out testing. If you
>>are
>> trying to do something like example shapefile_5.ncl, then what is
>>needed is
>> a coherent list of the coordinates of just the country boundary. In
>> shapefile_5.ncl, this was obtained through a shapefile that
>>contained *only*
>> that boundary.
>>
>> It would be easiest if you can find such a shapefile for your area
>>of
>> interest. If not however, you might be able to extract just the
>>required
>> polygon by subsetting the shapefile based upon a field in the
>>database.
>> Trying doing a "printVarSummary" on the shapefile file-variable.
>> Your
>> looking for a field that might encode boundary types, and most
>>importantly,
>> you want to make sure that the geometry_type attribute is "polygon".
>> Example shapefile_1.ncl shows how to query the shapefile's
>>database.
>>
>> If that's not clear and you want to send me your shapefiles, I'll
>>take a
>> look and see what might be possible.
>>
>> Good luck!
>> Rick
>>
>> On Dec 9, 2010, at 10:50 PM, shinn wong wrote:
>>
>> Dear Rick,
>> Thank you very much for your help. Now I can draw the country
>>outline, but
>> I've got another question. I want to check whether a point is within
>>the
>> country. So I use the function gc_inout(lat,lon,mrb_lat,mrb_lon).
>>But there
>> are points outside the country that show positive results. Does
>>anyone know
>> why? Thanks in advance for any given help.
>> Shinn
>> 2010/12/10 Rick Brownrigg <brownrig@ucar.edu>
>>
>>> Hi Shinn,
>>>
>>> The trouble stems from the single call to:
>>>
>>> line_data = gsn_add_polyline(wks, map_data, mrb_lon, mrb_lat,
>>>lnres)
>>>
>>> which will treat all the coordinates in mrb_lot/mrb_lat as a
>>>continuous
>>> line. You want to draw the line segments that make up the countries
>>> individually, as in shapefile_3.ncl. In particular, you want to
>>>replace the
>>> above call with something like:
>>>
>>> segments = f->segments
>>> geometry = f->geometry
>>> segsDims = dimsizes(segments)
>>> geomDims = dimsizes(geometry)
>>>
>>> geom_segIndex = f@geom_segIndex
>>> geom_numSegs = f@geom_numSegs
>>> segs_xyzIndex = f@segs_xyzIndex
>>> segs_numPnts = f@segs_numPnts
>>> numFeatures = geomDims(0)
>>>
>>> lines = new(segsDims(0),graphic) ; array to hold polylines
>>> segNum = 0 ; Counter for adding polylines
>>> do i=0, numFeatures-1
>>> startSegment = geometry(i, geom_segIndex)
>>> numSegments = geometry(i, geom_numSegs)
>>> do seg=startSegment, startSegment+numSegments-1
>>> startPT = segments(seg, segs_xyzIndex)
>>> endPT = startPT + segments(seg, segs_numPnts) - 1
>>> lines(segNum) = gsn_add_polyline(wks, map_data,
>>>mrb_lon(startPT:endPT), \
>>> mrb_lat(startPT:endPT),
>>>lnres)
>>> segNum = segNum + 1
>>> end do
>>> end do
>>>
>>>
>>> Check out the discussion at the top of the shapefiles page if the
>>>above
>>> code doesn't make sense. Hope that helps...
>>> Rick
>>>
>>>
>>>
>>>
>>>
>>> On Dec 9, 2010, at 2:00 AM, shinn wong wrote:
>>>
>>> Dear all,
>>> I've tried to plot the outline boundary of some countries using data
>>> obtained from DIVA-GIS. The code and the output file can be seen in
>>>the
>>> attachment. However, the output is quite different from that of
>>>example 4 in
>>> the shapefile example page. Does anyone know how can i just plot the
>>>outline
>>> of the country like those in example 4 instead of all the data
>>>points? How
>>> can i get the coordinate of the outline of the country? Many thanks.
>>> Shinn <aa1.ncl><shapefiles1.000001.png>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>>
>>
>>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Dec 13 10:52:14 2010

This archive was generated by hypermail 2.1.8 : Wed Dec 22 2010 - 16:10:23 MST