Re: about incompletely masked area with shapefile

This looks like a case of longitude values that go from 0 to 360 in your data, but from -180 to 180 in your shapefile outlines.

What do your "lon1d" values look like?

   print(min(lon1d))
   print(max(lon1d))

If you have any values greater than 180, then sometime before you assign "lon1d" to "data&lon", try adding:

   lon1d = where(lon1d.gt.180,lon1d-360,lon1d)

If this doesn't work, then it would really help if you could provide all the necessary data files. You can use our ftp:

http://www.ncl.ucar.edu/report_bug.shtml#HowToFTP

--Mary

On Jan 20, 2014, at 2:55 AM, dyjbean@gmail.com wrote:

> hi,
> i have met a question when i use shapefile file.
> <mask.png>
>
> in the middle panel, it's strange that the masked domain exists value, but in the bottom panel, the masked area have not any value.
> they have the same frame.
>
> my script used in terms of mask_12.ncl from www.ncl.ucar.edu/Applications/shapefiles.shtml , as follows:
>
> ;----------------------------------------------------------------------
> ; mask_12.ncl
> ;
> ; Concepts illustrated:
> ; - Using a worldwide shapefile to create a land/ocean mask
> ; - Masking a data array based on a geographical area
> ; - Attaching shapefile polylines to a map plot
> ; - Attaching lat/lon points to a map using gsn_coordinates
> ;----------------------------------------------------------------------
> ; Downloaded GSHHS shapefiles from:
> ;
> ; http://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/
> ;
> ; Used the "coarsest" one: "GSHHS_shp/c/GSHHS_c_L1.shp".
> ;----------------------------------------------------------------------
>
> 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"
>
> ;----------------------------------------------------------------------
> ; Function : create_mask_from_shapefile
> ;
> ; This function takes a shapefile of type "polygon" and a 2D data
> ; array that has coordinate arrays, and creates a mask array
> ; the same size as the data, that contains 0s and 1s, depending if
> ; the lat/lon values of the data are inside (1) or outside (0) of the
> ; shapefile polygons.
> ;----------------------------------------------------------------------
> undef("create_mask_from_shapefile")
> function create_mask_from_shapefile(data[*][*]:numeric,fname[1]:string)
> local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
> geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, ilat, ilon, i, j, \
> shp_lat, shp_lon, data_lat, data_lon, startSegment, numSegments, seg, \
> startPT, endPT, mask_val
> begin
> ;---Error checking
> if(.not.isfilepresent(fname)) then
> print("Error: create_mask_from_shapefile:")
> print(" '" + fname + "' doesn't exist.")
> print(" Mask array with all missing values will be returned.")
> return(new(dimsizes(data),integer,-999))
> end if
>
> ;---Check that "data" has 1D coordinate arrays.
> if(.not.isdimnamed(data,0).or..not.isdimnamed(data,1).or.\
> .not.iscoord(data,data!0).or..not.iscoord(data,data!1)) then
> print("Error: create_mask_from_shapefile:")
> print(" Input data doesn't have 1D coordinate arrays")
> print(" Mask array with all missing values will be returned.")
> return(new(dimsizes(data),integer,-999))
> end if
>
> ;---Open the shapefile
> f = addfile(fname,"r")
>
> ;---We can't use this routine to mask against point or line data
> if(f@geometry_type.ne."polygon") then
> print("Error: create_mask_from_shapefile: geometry_type attribute must be 'polygon'")
> print(" Mask array with all missing values will be returned.")
> return(new(dimsizes(data),integer,-999))
> end if
>
> ;---Read data off the shapefile
> segments = f->segments
> geometry = f->geometry
> segsDims = dimsizes(segments)
> geomDims = dimsizes(geometry)
>
> ;---Read global attributes
> geom_segIndex = f@geom_segIndex
> geom_numSegs = f@geom_numSegs
> segs_xyzIndex = f@segs_xyzIndex
> segs_numPnts = f@segs_numPnts
> numFeatures = geomDims(0)
>
> ;---Create mask array
> data_mask = new(dimsizes(data),integer,-999)
> mask_val = 1 ; 1's represent being inside a polygon
> data_mask = 0 ; 0's represent being outside a polygon
>
> ;---Read lat/lon values
> shp_lon = f->x
> shp_lat = f->y
> data_lat = data&$data!0$
> data_lon = data&$data!1$
> nlat = dimsizes(data_lat)
> nlon = dimsizes(data_lon)
>
> 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
> ;
> ; Messy code to make the nested loop go faster. Only check
> ; the area of data that is close to polygon that we are
> ; traversing. If you check data's full lat/lon array against
> ; every shapefile polygon, this loop will be extremely slow.
> ;
>
> ;
> ; Get approx. index values where data's lat/lon is close
> ; to shapefile's current polygon.
> ;
> iilt_beg = ind(data_lat.le.min(shp_lat(startPT:endPT)))
> iilt_end = ind(data_lat.ge.max(shp_lat(startPT:endPT)))
> iiln_beg = ind(data_lon.le.min(shp_lon(startPT:endPT)))
> iiln_end = ind(data_lon.ge.max(shp_lon(startPT:endPT)))
> ilt_beg = 0
> iln_beg = 0
> ilt_end = nlat-1
> iln_end = nlon-1
> if(.not.any(ismissing(iilt_beg))) then
> ilt_beg = iilt_beg(dimsizes(iilt_beg)-1)
> end if
> if(.not.any(ismissing(iilt_end))) then
> ilt_end = iilt_end(0)
> end if
> if(.not.any(ismissing(iiln_beg))) then
> iln_beg = iiln_beg(dimsizes(iiln_beg)-1)
> end if
> if(.not.any(ismissing(iiln_end))) then
> iln_end = iiln_end(0)
> end if
> ;
> ; Loop across subset of data's lat/lon and check each point
> ; to see if it is inside or outside of the shapefile polygon.
> ;
> do ilat = ilt_beg,ilt_end
> do ilon = iln_beg,iln_end
> if(data_mask(ilat,ilon).ne.mask_val.and.\
> gc_inout(data_lat(ilat),data_lon(ilon),\
> shp_lat(startPT:endPT),shp_lon(startPT:endPT))) then
> data_mask(ilat,ilon) = mask_val
> end if
> end do
> end do
> delete([/iilt_beg,iilt_end,iiln_beg,iiln_end/])
> end do
> end do
> return(data_mask)
> end
>
> ;----------------------------------------------------------------------
> ; Main code
> ;----------------------------------------------------------------------
> begin
>
> WRITE_MASK = True
> DEBUG = False
>
> dir="/public/home/douyj/wks_dyj/clm3.5/bld/clmrun_heihe_0.1d/"
> cdf_prefix="clmtest.clm2.h0.2008-08-21-00000"
> sfile=dir+cdf_prefix+".nc"
>
> ff=addfile(sfile,"r")
> ta=ff->TBOT(0,:,:)
> le1=ff->FGEV(0,:,:)
> le2=ff->FCEV(0,:,:)
> le3=ff->FCTR(0,:,:)
> lamda=(2.501-0.00236*(ta-273.15))*1000000
> data=(le1+le2+le3)/lamda*24*3600
> lat1d=ff->latixy(:,0)
> lon1d=ff->longxy(0,:)
> lat1d@units="degrees_north"
> lon1d@units="degrees_east"
>
> data!0="lat"
> data!1="lon"
> data&lat=lat1d
> data&lon=lon1d
> data@_FillValue=-9999
>
> shpfile="/public/home/douyj/wk/heihe_shapefile/heihe_hw.shp"
> heihe_mask=create_mask_from_shapefile(data,shpfile)
>
> ;---Mask "u" against land and ocean.
> u_land_mask = where(heihe_mask.eq.1,data,data@_FillValue)
> u_ocean_mask = where(heihe_mask.eq.0,data,data@_FillValue)
> copy_VarMeta(data,u_land_mask)
> copy_VarMeta(data,u_ocean_mask)
>
> ;---Start the graphics
> wks = gsn_open_wks("png","mask")
>
> res = True
>
> res@gsnMaximize = True ; maximize plot in frame
> res@gsnDraw = False ; don't draw plot yet
> res@gsnFrame = False ; don't advance frame yet
>
> res@cnFillOn = True
> res@cnLineLabelsOn = False
> res@cnLinesOn = False
>
> ;---Make sure both plots have same contour levels
> mnmxint = nice_mnmxintvl(min(data),max(data),25,False)
> res@cnLevelSelectionMode = "ManualLevels"
> res@cnMinLevelValF = mnmxint(0)
> res@cnMaxLevelValF = mnmxint(1)
> res@cnLevelSpacingF = mnmxint(2)
>
> res@lbLabelBarOn = False
> res@gsnAddCyclic = False
>
> res@mpFillOn = False
> res@mpOutlineOn = False
>
> res@gsnRightString = ""
> res@gsnLeftString = ""
>
> res@mpMinLatF = min(lat1d)
> res@mpMaxLatF = max(lat1d)
> res@mpMinLonF = min(lon1d)
> res@mpMaxLonF = max(lon1d)
> ;---Create plot of original data and attach shapefile outlines
> res@tiMainString = "Original data with shapefile outlines"
> map_data = gsn_csm_contour_map(wks,data,res)
> dum1 = gsn_add_shapefile_polylines(wks,map_data,shpfile,False)
>
> ;---Create plots of masked data
> res@tiMainString = "Original data masked against land"
> map_land_mask = gsn_csm_contour_map(wks,u_land_mask,res)
> res@tiMainString = "Original data masked against ocean"
> map_ocean_mask = gsn_csm_contour_map(wks,u_ocean_mask,res)
>
> if(DEBUG) then
> mkres = True
> ; mkres@gsMarkerSizeF = 0.007
> mkres@gsnCoordsAttach = True
> gsn_coordinates(wks,map_data,data,mkres)
> mkres@gsnCoordsNonMissingColor = "yellow"
> mkres@gsnCoordsMissingColor = "black"
> gsn_coordinates(wks,map_land_mask,u_land_mask,mkres)
> gsn_coordinates(wks,map_ocean_mask,u_ocean_mask,mkres)
> end if
>
> ;---Add shapefile outlines
> dum2 = gsn_add_shapefile_polylines(wks,map_land_mask,shpfile,False)
> dum3 = gsn_add_shapefile_polylines(wks,map_ocean_mask,shpfile,False)
>
> ;---Draw all three plots on one page
> pres = True
> pres@gsnMaximize = True
> pres@gsnPanelLabelBar = True
> gsn_panel(wks,(/map_data,map_land_mask,map_ocean_mask/),(/3,1/),pres)
>
> if(WRITE_MASK) then
> delete(ff) ; Close file before we open again.
> ;
> ; Make copy of file so we don't overwrite original.
> ; This is not necessary, but it's safer.
> ;
> new_cdf_file = cdf_prefix + "_with_mask.nc"
> system("/bin/cp " + sfile + " " + new_cdf_file)
> finout = addfile(new_cdf_file,"w")
> filevardef(finout, "land_mask", typeof(heihe_mask), (/ "lat", "lon" /) )
> finout->land_mask = (/heihe_mask/)
> end if
>
> end
>
> ;--------------------------------------------------------------------------
>
> my question is about how to get rid of the value in my masked area?
>
>
>
> thanks
>
>
> dyjbean@gmail.com
> _______________________________________________
> 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 Tue Jan 21 14:20:08 2014