[ncl-talk] Mask using shapefile
Mary Haley
haley at ucar.edu
Thu Mar 9 10:19:31 MST 2017
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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170309/b8546ea4/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mask_India_IND.000002.png
Type: image/png
Size: 96584 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170309/b8546ea4/attachment.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mask_India_IND.000001.png
Type: image/png
Size: 133931 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170309/b8546ea4/attachment-0001.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plot_mask.ncl
Type: application/octet-stream
Size: 6598 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170309/b8546ea4/attachment.obj
More information about the ncl-talk
mailing list