[ncl-talk] convert polyline shapefile to netcdf via ncl_convert2nc
Dave Allured - NOAA Affiliate
dave.allured at noaa.gov
Mon Apr 9 21:20:40 MDT 2018
Okay, I see now. A little more programming will be needed to count only
the unique SN values. There are several methods. Here is one way.
1. Before the loop, make sure the input list is sorted on "SN" so that all
identical SN values are grouped together. It doesn't matter whether it is
ascending or descending.
2. Before the loop, create two new grid arrays, in addition to count and
valavg. Let's call them inc (increment) and ucount. Set both to zero.
3. Inside the loop, on every valid cvalue, set inc(jl,il) = 1, for the
corresponding grid point only.
4. Inside the loop, after (3), whenever the next SN value is different, or
on the very last iteration, add the whole inc grid to the ucount grid.
Then set the inc grid back to all zeros. This has the effect of counting
each unique SN value only once in ucount.
At the end, ucount will contain the gridded unique counts that you want.
There may be an existing NCL example, but I am not aware of one.
--Dave
On Mon, Apr 9, 2018 at 8:03 PM, Lyndz <olagueralyndonmark429 at gmail.com>
wrote:
> Hi Sir Dennis and Sir Dave,
>
> Thank you for this. The code is actually doing what is was intended to do:
> counting the occurrence frequency of lat-lon points within a grid box.
>
> However, I only wanted to count the "unique" lat-lon points in each grid
> box. Each lat-lon points in the csv file has a corresponding "SN" value (SN
> column).
> If the lat-lon points, in each grid box, have the same "SN" value, then it
> should only be counted once.
>
> I'm not sure how to implement this in the code.
>
> Sincerely,
>
> Lyndz
>
>
> On Tue, Apr 10, 2018 at 10:58 AM, Dennis Shea <shea at ucar.edu> wrote:
>
>> I'll take Dave's example a bit further.
>>
>> Given that count, valavg and freq have coordinate variables. [Please
>> read.]
>>
>> Variable: count
>> Number of Dimensions: 2
>> Dimensions and sizes: [lat | 11] x [lon | 19]
>> Coordinates:
>> lat: [ 0..50]
>> lon: [90..180]
>> Number Of Attributes: 3
>> units :
>> long_name : Occurrence Count: vMax=35
>> _FillValue : 1e+20
>> ====================================================
>> test_lat = 20.0
>> test_lon = 130.0
>> print("---")
>> print("---> ("+test_lat+","+test_lon+")")
>> print("count ="+count({test_lat},{test_lon})) ; coordinate
>> subscripting {...}
>> print("valavg="+valavg({test_lat},{test_lon}))
>> print(" freq="+ freq({test_lat},{test_lon}))
>> print("---")
>>
>> (0) ---> (20,130)
>> (0) count =86
>> (0) valavg=76.1047
>> (0) freq=7.83242
>>
>> =========================================================
>> See attached
>>
>> On Mon, Apr 9, 2018 at 7:12 PM, Dave Allured - NOAA Affiliate <
>> dave.allured at noaa.gov> wrote:
>>
>>> Oops, there is a bug in my print statement. Change "lines" to "data":
>>>
>>> print (count(jl,il) + " " + data(n))
>>>
>>>
>>> On Mon, Apr 9, 2018 at 6:22 PM, Dave Allured - NOAA Affiliate <
>>> dave.allured at noaa.gov> wrote:
>>>
>>>> Lyndz,
>>>>
>>>> Try this to find the problem with your counting method. Below this
>>>> line:
>>>>
>>>> valavg(jl,il) = valavg(jl,il) + cval(n)
>>>>
>>>> Add these code lines:
>>>>
>>>> test_lat = 20.0
>>>> test_lon = 130.0
>>>> jtest = toint((test_lat-latS)/dlat)
>>>> itest = toint((test_lon-lonW)/dlon)
>>>>
>>>> if (il .eq. itest .and. jl .eq. jtest) then
>>>> print (count(jl,il) + " " + lines(n))
>>>> end if
>>>>
>>>> In the code above, change test_lat and test_lon to the coordinates of
>>>> any single grid box that you want to check. Then run the program. This
>>>> will print every input record that gets counted for that one grid box.
>>>>
>>>> For example, I get 86 counts for the grid box example above, i..e. 20N,
>>>> 130E.
>>>>
>>>> --Dave
>>>>
>>>>
>>>> On Sun, Apr 8, 2018 at 9:04 PM, Lyndz <olagueralyndonmark429 at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi Sir Dennis,
>>>>>
>>>>> I tried your suggestion. I was able to run the script but I'm getting
>>>>> a very large count. When I counted it manually, the count does not exceed
>>>>> 10. But I'm getting around 88.
>>>>> Attached are the output image, csv input file, and modified script.
>>>>>
>>>>> I only wanted to count the TC with Vmax greater than 35knts. So I
>>>>> tried masking first the data.
>>>>>
>>>>> I might be missing one important step here.
>>>>>
>>>>> Many thanks for the help.
>>>>>
>>>>> Sincerely,
>>>>>
>>>>> Lyndz
>>>>>
>>>>>
>>>>> On Mon, Apr 9, 2018 at 3:17 AM, Dennis Shea <shea at ucar.edu> wrote:
>>>>>
>>>>>> http://www.ncl.ucar.edu/Applications/binning.shtml
>>>>>>
>>>>>> See Example 2:
>>>>>>
>>>>>> freq = count
>>>>>> freq = (count/npts)*100 ; comment this
>>>>>>
>>>>>>
>>>>>> On Sat, Apr 7, 2018 at 9:05 PM, Lyndz <olagueralyndonmark429 at gmail.c
>>>>>> om> wrote:
>>>>>>
>>>>>>> Hi Sir Rick,
>>>>>>>
>>>>>>> Thank you for the fast response.
>>>>>>> I wanted to create a 5degree by 5degree grid NetCDF file from the
>>>>>>> shapefile, where each grid box contains the "count"(frequency) of unique TC
>>>>>>> tracks per category.
>>>>>>> I'm just wondering if this is possible to do in NCL.
>>>>>>>
>>>>>>> Sincerely,
>>>>>>>
>>>>>>> Lyndz
>>>>>>>
>>>>>>>
>>>>>>> On Sun, Apr 8, 2018 at 3:34 AM, Rick Brownrigg <brownrig at ucar.edu>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hi Lyndz,
>>>>>>>>
>>>>>>>> I'm not sure there is a "correct" way to do this. NetCDF is really
>>>>>>>> good at representing regular arrays of data, whereas polylines/polygons
>>>>>>>> tend to have varying numbers of coordinate pairs -- what would lat and lon
>>>>>>>> variables look like in that case, and what would be the meaning of lat/lon
>>>>>>>> dimensions? NCL makes a shapefile *look* like a NetCDF file by packing the
>>>>>>>> coordinate information for all features into the x/y variables, but then
>>>>>>>> one has to utilize the geometry and segments variables to unpack
>>>>>>>> coordinates for each feature. The shapefiles examples page show many
>>>>>>>> examples of doing this:
>>>>>>>>
>>>>>>>> http://ncl.ucar.edu/Applications/shapefiles.shtml
>>>>>>>>
>>>>>>>> If you are working with NCL, you are probably better off leaving
>>>>>>>> your data as a shapefile*. *I take it however that you have other
>>>>>>>> reasons for wanting a NetCDF file?
>>>>>>>>
>>>>>>>> Finally, I might comment that the conversion result may not be
>>>>>>>> correct -- its suspect to me that the "num_points" value is exactly twice
>>>>>>>> the value of "num_features" -- I wonder if maybe just the end-points of the
>>>>>>>> tracks where what got captured in the conversion?
>>>>>>>>
>>>>>>>> Rick
>>>>>>>>
>>>>>>>> On Fri, Apr 6, 2018 at 11:30 PM, Lyndz <
>>>>>>>> olagueralyndonmark429 at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> Dear NCL experts,
>>>>>>>>>
>>>>>>>>> I would like to convert the following shapefile to a netcdf file.
>>>>>>>>> I created the shapefile from a csv file containing jtwc tc tracks
>>>>>>>>> (see attached csv2shp.py)
>>>>>>>>> Also attached is the csv file.
>>>>>>>>>
>>>>>>>>> When I used the ncl_convert2nc the netcdf file has no lat-lon
>>>>>>>>> dimension. Here's the output of the ncdump. Is it also possible to divide
>>>>>>>>> the TC categories by dividing the Vmax values?
>>>>>>>>>
>>>>>>>>> I'll appreciate any suggestion on how to do this correctly.
>>>>>>>>>
>>>>>>>>> netcdf par_jtwc_above_ts_1979-1993 {
>>>>>>>>> dimensions:
>>>>>>>>> geometry = 2 ;
>>>>>>>>> segments = 2 ;
>>>>>>>>> num_features = 1050 ;
>>>>>>>>> num_segments = 1050 ;
>>>>>>>>> num_points = 2100 ;
>>>>>>>>> variables:
>>>>>>>>> int geometry(num_features, geometry) ;
>>>>>>>>> int segments(num_segments, segments) ;
>>>>>>>>> double x(num_points) ;
>>>>>>>>> double y(num_points) ;
>>>>>>>>> int SN(num_features) ;
>>>>>>>>> int CY(num_features) ;
>>>>>>>>> int Y1(num_features) ;
>>>>>>>>> int M1(num_features) ;
>>>>>>>>> int D1(num_features) ;
>>>>>>>>> int H1(num_features) ;
>>>>>>>>> int VMax1(num_features) ;
>>>>>>>>> int Y2(num_features) ;
>>>>>>>>> int M2(num_features) ;
>>>>>>>>> int D2(num_features) ;
>>>>>>>>> int H2(num_features) ;
>>>>>>>>> int VMax2(num_features) ;
>>>>>>>>>
>>>>>>>>> // global attributes:
>>>>>>>>> :segs_numPnts = 1 ;
>>>>>>>>> :segs_xyzIndex = 0 ;
>>>>>>>>> :geom_numSegs = 1 ;
>>>>>>>>> :geom_segIndex = 0 ;
>>>>>>>>> :geometry_type = "polyline" ;
>>>>>>>>> :layer_name = "par_jtwc_above_ts_1979-1993" ;
>>>>>>>>> :creation_date = "Sat Apr 7 14:11:24 JST 2018" ;
>>>>>>>>> :NCL_Version = "6.4.0" ;
>>>>>>>>> :system = "Linux localhost.localdomain 3.10.0-327.36.3.el7.x86_64
>>>>>>>>> #1 SMP Mon Oct 24 16:09:20 UTC 2016 x86_64 x86_64 x86_64 GNU/Linux" ;
>>>>>>>>> :Conventions = "None" ;
>>>>>>>>> :title = "NCL: convert-OGR-to-netCDF" ;
>>>>>>>>> }
>>>>>>>>>
>>>>>>>>> Sincerely,
>>>>>>>>>
>>>>>>>>> Lyndz
>>>>>>>>>
>>>>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180409/7f0f4fc3/attachment.html>
More information about the ncl-talk
mailing list