[ncl-talk] Mask using shapefile
Ipshita Majhi
ipmajhi at alaska.edu
Thu Mar 9 11:52:08 MST 2017
Thank you so much for this mask and the plot.
On Thu, Mar 9, 2017 at 8:19 AM, Mary Haley <haley at ucar.edu> wrote:
> Ipshita,
>
> Thanks for the file.
>
> There were a few errors in your script.
>
> The big one is that you don't have an "end" statement for your function.
> You instead put two "end" statements at the very end of the script.
>
> You need to remove one of the "end" statements at the end of the script
> and add one right after the line:
>
> return(data_mask)
>
> I think this might get your code to work again.
>
> I was able to run your script with this fix, but I first had to correct
> the filename and the variable name, because you gave me a different file
> than what is used in the script.
>
> Also, you have:
>
> res at tiMainString = filename
>
> but "filename" is not defined. I defined "filename" and changed your
> "addfile" line to this:
>
> filename = "b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECSL.070001-079912.nc"
> f = addfile(filename,"r")
>
>
> I don't have the India_adm files, so I changed this code to use IND_adm
> instead, which I already had.
>
>
> --Mary
>
> On Wed, Mar 8, 2017 at 5:44 PM, Ipshita Majhi <ipmajhi at alaska.edu> wrote:
>
>> Hi Mary,
>>
>> Here are the errors I am getting after incorporating what you
>> recommended.
>>
>> (0) gsn_csm_map_ce: Warning: you set mpMaxLonF to a value > 180, but
>> (0) didn't set mpCenterLonF. Setting mpCenterLonF to
>> 179.375
>> ERROR 10: Pointer 'hFeat' is NULL in 'OGR_F_GetGeometryRef'.
>>
>> ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryCount'.
>>
>> ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetPointCount'.
>>
>> fatal:["Execute.c":5734]:variable (geometry) is not in file (f)
>> fatal:["Execute.c":7743]:Execute: Error occurred at or near line 26 in
>> file shapefile_India.ncl
>>
>> fatal:["Execute.c":7743]:Execute: Error occurred at or near line 145 in
>> file shapefile_India.ncl
>> **************************************
>> Is it because my version of NCL is 6.1.0?
>>
>> Here is my file that you had requested
>> b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECSL.07000...
>> <https://drive.google.com/a/alaska.edu/file/d/0BxwHbkyQt46bU2J0VGwyREh0eTA/view?usp=drive_web>
>>
>>
>> On Wed, Mar 8, 2017 at 2:36 PM, Mary Haley <haley at ucar.edu> wrote:
>>
>>> Hi,
>>>
>>> Can you provide the b.e11 NetCDF file? You can email me offline about
>>> this if you want.
>>>
>>> I do recognize the code you included because it's code that I wrote, but
>>> it's hard to debug without the data file.
>>>
>>> Also, did you look at the "shapefile_mask_data" function on the
>>> shapefiles example page?
>>>
>>> In particular, you could look at example shapefile_11.ncl:
>>>
>>> http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex11
>>>
>>> You can first try it with minimal options to see what happens:
>>>
>>> opt = True
>>> opt at debug = True
>>> ts_mask = shapefile_mask_data(ts,India_shp,opt)
>>> printVarSummary(ts_mask)
>>> printMinMax(ts_mask,0)
>>>
>>> I can't remember, but you might need to call this too:
>>>
>>> copy_VarMeta(ts,ts_mask) ; Copy all metadata
>>>
>>>
>>> --Mary
>>>
>>>
>>>
>>>
>>>
>>> On Wed, Mar 8, 2017 at 12:59 PM, Ipshita Majhi <ipmajhi at alaska.edu>
>>> wrote:
>>>
>>>> Dear NCL Users,
>>>>
>>>> I am trying to create a mask for India, using shapefile. I am not
>>>> getting any error but there is no plot begin created. I will be grateful if
>>>> someone could guide me about it. Here is the code:
>>>>
>>>> *********************************************************
>>>> ;************************************************
>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
>>>> ;************************************************
>>>> ;Read data to plot and mask
>>>> ;************************************************
>>>> function mask_data_with_India_country(country_code,data)
>>>> begin
>>>>
>>>> mask_start_time = get_cpu_time() ; For timing results
>>>>
>>>> ;---Convert rectilinear grid to 2D grid, but laid out as 1D array.
>>>> dims = dimsizes(data)
>>>> lat1d = ndtooned(conform_dims(dims,data&$data!0$,0))
>>>> lon1d = ndtooned(conform_dims(dims,data&$data!1$,1))
>>>>
>>>> shapefile_name = "IND_adm0.shp"
>>>> ;---Open shapefile and read lat/lon values.
>>>> sfilename = "IND_adm0.shp"
>>>> f = addfile(sfilename,"r")
>>>> segments = f->segments
>>>> geometry = f->geometry
>>>> segsDims = dimsizes(segments)
>>>> geomDims = dimsizes(geometry)
>>>>
>>>> ;---Read global attributes
>>>> geom_segIndex = f at geom_segIndex
>>>> geom_numSegs = f at geom_numSegs
>>>> segs_xyzIndex = f at segs_xyzIndex
>>>> segs_numPnts = f at segs_numPnts
>>>> numFeatures = geomDims(0)
>>>>
>>>> lon = f->x
>>>> lat = f->y
>>>> nlatlon = dimsizes(lat)
>>>>
>>>> min_lat = min(lat)
>>>> max_lat = max(lat)
>>>> min_lon = min(lon)
>>>> max_lon = max(lon)
>>>>
>>>> print("==================================================")
>>>> print("Shapefile : " + sfilename)
>>>> print("min_lat " + min_lat)
>>>> print("max_lat " + max_lat)
>>>> print("min_lon " + min_lon)
>>>> print("max_lon " + max_lon)
>>>>
>>>> ii_latlon = ind(min_lat.le.lat1d.and.lat1d.le.max_lat.and.\
>>>> min_lon.le.lon1d.and.lon1d.le.max_lon)
>>>> nii = dimsizes(ii_latlon)
>>>> print(nii + " values to check")
>>>> print(numFeatures + " feature(s) to traverse with a maximum of " + \
>>>> nlatlon + " points")
>>>>
>>>>
>>>> ;---Create array to hold new data mask, and set all values to 0
>>>> initially.
>>>> data_mask_1d = new(dimsizes(lat1d),integer)
>>>> data_mask_1d = 0
>>>>
>>>> ;
>>>> ; This is the loop that checks every point in lat1d/lon1d to see if it
>>>> ; is inside or outside of the country. If it is inside, then
>>>> data_mask_1d
>>>> ; will be set to 1.
>>>> ;
>>>> ikeep = 0 ; Counter to see how many points were found inside the
>>>> country
>>>> do n=0,nii-1
>>>> ii = ii_latlon(n)
>>>> is_inside = False
>>>> i = 0
>>>> do while(.not.is_inside.and.i.lt.numFeatures)
>>>> 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
>>>> if(data_mask_1d(ii).ne.1.and.\
>>>> gc_inout(lat1d(ii),lon1d(ii),\
>>>> lat(startPT:endPT),lon(startPT:endPT))) then
>>>> data_mask_1d(ii) = 1
>>>> ikeep = ikeep+1
>>>> is_inside = True
>>>> continue
>>>> end if
>>>> end do
>>>> i = i + 1
>>>> end do
>>>> end do
>>>> print(ikeep + " values kept")
>>>> print("==================================================")
>>>>
>>>> ; Create a 2D data array of same size as original data,
>>>> ; but with appropriate values masked.
>>>> ;
>>>> data_mask = (where(onedtond(data_mask_1d,dims).eq.1,data,\
>>>> data at _FillValue))
>>>> copy_VarMeta(data,data_mask) ; Copy all metadata
>>>>
>>>> ;---Print timings
>>>> mask_end_time = get_cpu_time()
>>>> print("Elapsed time in CPU second for 'mask_data_with_India_country'
>>>> = " + \
>>>> (mask_end_time-mask_start_time))
>>>>
>>>> return(data_mask)
>>>>
>>>> ;----------------------------------------------------------------------
>>>> ; Main code
>>>> ;----------------------------------------------------------------------
>>>> begin
>>>> ;---Shapefile to use for masking. Using "Bolivia" here (BOL).
>>>> India_code = "India"
>>>> India_name = India_code + "_adm"
>>>> India_shp = India_name + "/" + India_name + "0.shp"
>>>>
>>>> ;---Read lat/lon off shapefile
>>>> s = addfile(India_shp,"r")
>>>> slat = s->y
>>>> slon = s->x
>>>>
>>>> ;---Read precipitation data to contour and mask
>>>> f = addfile("b.e11.B1850C5CN.f09_
>>>> g16.005.cam.h0.PRECC.070001-079912.nc","r")
>>>> ts = f->PRECC(0,:,:)
>>>> printVarSummary(ts)
>>>>
>>>> ;---Start the graphics
>>>> wtype = "x11" ; "png"
>>>> ; wtype at wkWidth = 2500 ; use for "png" or "x11"
>>>> ; wtype at wkHeight = 2500
>>>>
>>>> wks = gsn_open_wks(wtype,"mask_India_"+India_code)
>>>>
>>>> res = True ; plot mods desired
>>>>
>>>> res at cnFillOn = True ; turn on color
>>>> res at cnFillMode = "RasterFill"
>>>> res at cnLinesOn = False ; turn off lines
>>>> res at cnLineLabelsOn = False ; turn off labels
>>>>
>>>> res at mpMinLatF = min(ts&lat)
>>>> res at mpMaxLatF = max(ts&lat)
>>>> res at mpMinLonF = min(ts&lon)
>>>> res at mpMaxLonF = max(ts&lon)
>>>>
>>>> ;---Create global plot of original data
>>>> res at tiMainString = filename
>>>> plot = gsn_csm_contour_map(wks,ts, res)
>>>>
>>>> ;---Mask "ts" against India shapefile outlines
>>>> ts_mask = mask_data_with_India_country(India_code,ts)
>>>> printVarSummary(ts_mask)
>>>>
>>>> ;---Set some additional resources for the second set of plots
>>>>
>>>> res at gsnDraw = False ; don't draw yet
>>>> res at gsnFrame = False ; don't advance frame
>>>> yet
>>>>
>>>> ;---Pick "nice" contour levels for both plots
>>>> mnmxint = nice_mnmxintvl( min(ts_mask), max(ts_mask), 18, False)
>>>> res at cnLevelSelectionMode = "ManualLevels"
>>>> res at cnMinLevelValF = mnmxint(0)
>>>> res at cnMaxLevelValF = mnmxint(1)
>>>> res at cnLevelSpacingF = mnmxint(2)/2.
>>>>
>>>> res at mpMinLatF = min(slat)
>>>> res at mpMaxLatF = max(slat)
>>>> res at mpMinLonF = min(slon)
>>>> res at mpMaxLonF = max(slon)
>>>>
>>>> res at gsnRightString = ""
>>>> res at gsnLeftString = ""
>>>>
>>>> res at lbLabelBarOn = False
>>>>
>>>> ;---Create plot of original data
>>>> res at tiMainString = "original data zoomed in"
>>>> plot_orig = gsn_csm_contour_map(wks,ts, res)
>>>>
>>>> ;---Create plot of masked data
>>>> res at tiMainString = "with country mask"
>>>> plot_mask = gsn_csm_contour_map(wks,ts_mask, res)
>>>>
>>>> ;---Add shapefile outlines to both plots
>>>> lnres = True
>>>> lnres at gsLineColor = "NavyBlue"
>>>>
>>>> id_orig = gsn_add_shapefile_polylines(wks,plot_orig,India_shp,lnres)
>>>> id_mask = gsn_add_shapefile_polylines(wks,plot_mask,India_shp,lnres)
>>>>
>>>> ;---Panel both plots on one page
>>>> pres = True
>>>> pres at txString = ts at long_name + " (" + ts at units + ")"
>>>> pres at gsnMaximize = True
>>>> pres at gsnPanelLabelBar = True
>>>> gsn_panel(wks,(/plot_orig,plot_mask/),(/1,2/),pres)
>>>> end
>>>> end
>>>>
>>>>
>>>> --
>>>> ************************************************************
>>>> *******************************************************
>>>> "I slept and dreamt that life was joy. I awoke and saw that life was
>>>> service. I acted and behold, service was joy." - Rabindranath Tagore
>>>> ************************************************************
>>>> ********************************************************
>>>>
>>>> Ipshita Majhi
>>>> PhD Candidate
>>>> University of Alaska , Fairbanks
>>>> Atmospheric Science Department
>>>> (907)978-4220 <(907)%20978-4220> ipmajhi at alaska.edu
>>>>
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> ncl-talk at ucar.edu
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>>
>>>
>>
>>
>> --
>> ************************************************************
>> *******************************************************
>> "I slept and dreamt that life was joy. I awoke and saw that life was
>> service. I acted and behold, service was joy." - Rabindranath Tagore
>> ************************************************************
>> ********************************************************
>>
>> Ipshita Majhi
>> PhD Candidate
>> University of Alaska , Fairbanks
>> Atmospheric Science Department
>> (907)978-4220 <(907)%20978-4220> ipmajhi at alaska.edu
>>
>
>
--
*******************************************************************************************************************
"I slept and dreamt that life was joy. I awoke and saw that life was
service. I acted and behold, service was joy." - Rabindranath Tagore
********************************************************************************************************************
Ipshita Majhi
PhD Candidate
University of Alaska , Fairbanks
Atmospheric Science Department
(907)978-4220 ipmajhi at alaska.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170309/bcba596e/attachment-0001.html
More information about the ncl-talk
mailing list