[ncl-talk] Using Shapefiles to create subset from global forecast datasets

Klemm, Toni toni at ou.edu
Wed May 24 15:19:34 MDT 2017


Hi Rick,

Thanks for getting back to me!

I got it to work in the meantime using masking example 9 (https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl), which looks similar to example 4 on the page you linked (http://ncl.ucar.edu/Applications/Scripts/shapefiles_4.ncl). The masking worked, but for some reason there is large wedge with missing data right in the middle of my plot now. I imagine it has somewhere to do with the “if (MASK_INSIDE) then” loop in my script below, but I can’t figure it out. Do you have any suggestions? I also attached my original data file, the resulting datafile (_shp.nc) and plot, and the .nc file of my study area. This shapefile was created in ArcMap 10.3 from the states.* files in the NCL shapefile set<https://www.ncl.ucar.edu/Applications/Data/#shp> and converted to .nc using ncl_convert2nc.

Again, thank you very much!

Best,
Toni

Script:

; *****************************************************************
; Create map of GFDL-OBS bias averages by FORECAST MONTH 1980-2014
;
; Use SHAPEFILE OF STUDY REGION to define subset and plot
; *****************************************************************

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"

; ************************************************************

begin

do target_month = 1,12  ; do-loop through the 12 calendar months

  if (target_month .lt. 10) then
    target_month_zero = "0" + target_month ; turns "1"..."9" into "01"..."09" to match filename structure (months have 2 digits"
  end if

 if (target_month .gt. 9) then
    target_month_zero = target_month ; leaves "10..."12" unchanged
 end if

 print("****  START: target_month " + target_month + "  ****")

  ; *****************************************************
  ; READ IN FILE INFORMATION BY TARGET MONTH
  ; *****************************************************

; *****
print("1 - Read in file information")

  bias_f = addfile("~/GFDL/TAS/DATAFILES/monthly/5b-4alt_GFDL_TAS_bias/avg_by_target_month/stack/GFDL_TAS_bias_1980-2014_target_" + target_month_zero + ".nc","r")
  print("Finished: 1/8 Read in multiple files")


  ; ***********************************************
  ; READ IN THE VARIABLE (precipitation)
  ; ***********************************************

  print("2 - Read in variable and dimensions and average all years")
  lat = (/bias_f->lat/)
  lon = (/bias_f->lon/)

  printMinMax(lat, False)
  printMinMax(lon, False)

  nlat = dimsizes(lat)
  lat!0 = "lat"
  lat&lat = lat
  lat at long_name = "latitude"
  lat at units = "degrees_north"

  mlon = dimsizes(lon)
  lon!0 = "lon"
  lon&lon = lon
  lon at long_name = "longitude"
  lon at units = "degrees_east"

  GFDL_bias_3D = rm_single_dims(bias_f->GFDL_OBS_bias) ; extract precipitation into new variable, output is mm/day
  GFDL_bias_2D = dim_avg_n(GFDL_bias_3D,0)   ; averaging all bias for GFDL and LTM made in that month for 1985-2010

  min_bias = decimalPlaces(min(GFDL_bias_2D),2,True)
  max_bias = decimalPlaces(max(GFDL_bias_2D),2,True)
  avg_bias = decimalPlaces(avg(GFDL_bias_2D),2,True)

;  printVarSummary(GFDL_bias_2D)

;  GFDL_bias_2D at _FillValue = 1.0E20

;  printVarSummary(GFDL_bias_2D)

;  printVarSummary(GFDL_bias_3D)
;  printVarSummary(GFDL_bias_2D)


; *****************************************
; ***** Read in shapefile information ***** (source: https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl)

  MASK_INSIDE = True    ; Whether to mask data inside or outside the given geographical area.

  shp = addfile("~/States_shapefile_data/NCL/4_states_merged.shp","r")
  shp_lat = shp->y
  shp_lon = shp->x+360          ; adding 360 to get shapefile longitudes in same value range as data longitudes
;  printVarSummary(shp_lat)
;  printVarSummary(shp_lon)
  printMinMax(shp_lat, False)
  printMinMax(shp_lon, False)

  nshp = dimsizes(shp_lon)

  min_shp_lat = min(shp_lat)
  max_shp_lat = max(shp_lat)
  min_shp_lon = min(shp_lon)
  max_shp_lon = max(shp_lon)


; ***** Create new subset based on shapefile ***** (source: https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl)

  if(MASK_INSIDE) then
;---Start with data filled in.
    data_mask = GFDL_bias_2D
  else
;---Start with data all missing
  data_mask = new(dimsizes(GFDL_bias_2D),typeof(GFDL_bias_2D),GFDL_bias_2D at _FillValue)
  copy_VarCoords(GFDL_bias_2D,data_mask)
  end if

; ***** Determine the boundary lat/lons of the subset region within the lat/lon box
  ilt_mn = ind(min_shp_lat.gt.lat)
  ilt_mx = ind(max_shp_lat.lt.lat)
  iln_mn = ind(min_shp_lon.gt.lon)
  iln_mx = ind(max_shp_lon.lt.lon)
  ilt1   = ilt_mn(dimsizes(ilt_mn)-1)    ; Start of lat box
  iln1   = iln_mn(dimsizes(iln_mn)-1)    ; Start of lon box
  ilt2   = ilt_mx(0)                     ; End of lat box
  iln2   = iln_mx(0)                     ; End of lon box


  if (MASK_INSIDE) then
; ***** Put missing values in the areas that we want masked. (source: https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl)
    do ilt=ilt1,ilt2
      do iln=iln1,iln2
        if(gc_inout(lat(ilt),lon(iln),shp_lat,shp_lon)) then
          data_mask(ilt,iln) = data_mask at _FillValue                 ; replace data inside and outside study region with missing values
        end if
      end do
    end do
  else
; ***** Put data back in the areas that we don't want masked. (source: https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl)
    do ilt=ilt1,ilt2
      do iln=iln1,iln2
        if(gc_inout(lat(ilt),lon(iln),shp_lat,shp_lon)) then
          data_mask(ilt,iln) = GFDL_bias_2D(ilt,iln)                ; replace missing values *inside* study region with data again
        end if
      end do
    end do
  end if


  ; ***********************************************
  ; WRITE THE AVERAGES INTO NEW netCDF FILE
  ; ***********************************************

  print("3 - Write the averages into new netCDF file")

    netCDF1 = False ; Output format is NetCDF

  if (netCDF1) then
    diro = "~/GFDL/TAS/DATAFILES/monthly/5b-4alt_GFDL_TAS_bias/avg_by_target_month/nc_shp/" ; output in subdirectory "subset" within input directory

    filename = "GFDL_TAS_bias_1980-2014_target_month_" + target_month_zero + "_shp" ; output filename for subset using initialization year and month to differentiate
    filo = filename +".nc"

  end if

  setfileoption("nc","Format","LargeFile")

  if (netCDF1) then
    system("/bin/rm -f "+ diro + filo) ; removes old files with the same name if they are present
    ncdf = addfile(diro + filo,"c")

    setfileoption(ncdf,"DefineMode",True)

    ; create attributes
    fAtt = True
    fAtt at title = "Data: GFDL FLOR-B01 for the Contiguous U.S., Forecast Month " + target_month_zero
    fAtt at source = "Vecchi et al. 2014"
    fAtt at Conventions = "None"
    fAtt at creation_date = systemfunc("date")

    fileattdef(ncdf,fAtt)

    dimNames = (/"target_month","lat","lon"/)
    dimSizes = (/dimsizes(target_month),dimsizes(lat),dimsizes(lon)/)
    dimUnlim = (/True,False,False/)

    filedimdef(ncdf,dimNames,dimSizes,dimUnlim)

    filevardef(ncdf,"target_month",typeof(target_month),(/"target_month"/))
    filevardef(ncdf,"lat",typeof(lat),(/"lat"/))
    filevardef(ncdf,"lon",typeof(lon),(/"lon"/))
    filevardef(ncdf,"data_mask",typeof(data_mask),(/"lat","lon"/))

    ; now write all the variables to the file

    ncdf->target_month = (/target_month/)
    ncdf->lat = (/lat/)
    ncdf->lon = (/lon/)
    ncdf->data_mask = (/data_mask/)

  end if


  ; ***********************************************
  ; CREATE MAP
  ; ***********************************************

  ; Define boundaries for map
  data_mask!0 = "lat" ; latitude information come from variable lat_subdomain
  data_mask&lat = lat
  data_mask!1 = "lon" ; longitude information come from variable lon_subdomain
  data_mask&lon = lon
  data_mask at long_name = "Range: " + min_bias + " - " + max_bias + " mm/day"    ; map subtitle, printed on the top left
  data_mask at units = "Avg.: " + avg_bias + " mm/day" ; map units, printed on the top right

  diro = "~/GFDL/TAS/DATAFILES/monthly/5b-4alt_GFDL_TAS_bias/avg_by_target_month/png_shp/" ; output directory
  wks = gsn_open_wks("png", diro + "GFDL_TAS_bias_1980-2014_target_month_" + target_month_zero + "_shp") ; map name and file type

  cmap = read_colormap_file("BlueGreen14")

  res                      = True ; plot mods desired
  res at tiMainString         = "Avg. GFDL bias, Forecast Month " + target_month_zero + ", 1980-2014, no bias correction" ; main title
  res at cnFillOn             = True ; turn on color fill, works with gsn_csm_map() function 3 lines down
  res at gsnAddCyclic         = False
  res at cnLinesOn            = False ; contour lines of the map fill on or off
  res at mpProjection         = "Mercator" ; Mercator map projection
  res at cnFillPalette        = cmap ; reverse color map

  res at mpLimitMode    = "LatLon" ; define plotted area by lats and lons, namely those of the selection (20-50N, 230-300E)
  res at mpMaxLatF      = 41.5
  res at mpMinLatF      = 25
  res at mpMaxLonF      = 267
  res at mpMinLonF      = 250

  res at mpOutlineBoundarySets = "GeophysicalAndUSStates" ; lines and boundaries of USA and US States
  res at mpGeophysicalLineThicknessF = 2             ; thickness of the USA lines
  res at mpUSStateLineThicknessF     = 2               ; thickness of the state lines

  res at mpSpecifiedFillColors  = (/0,100/)             ; fill with background color, works with gsn_map() 2 lines down

;  res at cnLevelSelectionMode = "ManualLevels"              ; manually set the contour levels with the following 3 resources
;  res at cnMinLevelValF  = -2                             ; set the minimum contour level
;  res at cnMaxLevelValF  = 4                           ; set the maximum contour level
;  res at cnLevelSpacingF = 0.25                             ; set the interval between contours


  ; ***********************************************
  ; PLOT MAP
  ; ***********************************************

  print ("5 - Plot map")

  plot = gsn_csm_contour_map(wks,data_mask, res) ; universal map projection, define in res at mpProjection

  print("**** ALL DONE: target_month " + target_month + " ****")

  delete(bias_f)    ; delete before the next loop or program will fail
  delete(GFDL_bias_3D) ; delete before the next loop or program will fail
  delete(GFDL_bias_2D) ; delete before the next loop or program will fail
  delete(data_mask)   ; delete before the next loop or program will fail

end do ; end the loop through the 12 calendar months (1-12)

end

Toni Klemm
Ph.D. Candidate
The University of Oklahoma
South Central Climate Science Center
phone: 405 325 1272
www.toni-klemm.de<http://www.toni-klemm.de>
Early Career Climate Forum<https://www.eccforum.org>
Twitter: @toniklemm<https://twitter.com/ToniKlemm>


[cid:848bb57a-8d3f-4131-8b70-e31ed6753761 at namprd03.prod.outlook.com]



On May 24, 2017, at 12:35 PM, Rick Brownrigg <brownrig at ucar.edu<mailto:brownrig at ucar.edu>> wrote:

Hi Toni,

Sorry for the delay in response. There are several examples of using shapefiles to mask out data, if that's what are wanting to do. In particular, look at examples 4, 5, and 21 at:

http://ncl.ucar.edu/Applications/shapefiles.shtml<https://urldefense.proofpoint.com/v2/url?u=http-3A__ncl.ucar.edu_Applications_shapefiles.shtml&d=DwMFaQ&c=qKdtBuuu6dQK9MsRUVJ2DPXW6oayO8fu4TfEHS8sGNk&r=trcSnQNnUT2ddBXcpd358A&m=_2unS1RvDS7egGpkWknJW0V9KMHhv2GWDckbCG7ZOcQ&s=VqnvWQCKomt1QM_SAC2MNal5QUSSlnZmBwXdI3_KjO0&e=>

II think there's a utility script they make use that you'll need to download.  Hopefully this gets you at least started -- please post back the list if you have other questions.

Rick

On Mon, May 22, 2017 at 10:10 PM, Klemm, Toni <toni at ou.edu<mailto:toni at ou.edu>> wrote:
Good evening,

I am analyzing seasonal climate forecasts skill for the U.S. Great Plains, specifically Texas, Oklahoma, Kansas, and Colorado. So far I’ve been using lat/lon to define my study region. However, this rectangular box includes parts of Mexico and the Gulf of Mexico that are outside my study region and which have extreme outliers that distort my results.

To get more accurate results, I’d like to use a shapefile of my 4 states (TX_OK_KS_CO.shp, created in QGIS 2.18 from a U.S. states shapefile) to define my study region. However, I’m not sure how to interpret and use the results I’m getting. I used   “ncl_convert2nc TX_OK_KS_CO.shp -nc2c“  to create   “TX_OK_KS_CO.nc“ .

“printVarSummary(TX_OK_KS_CO.nc)” and ”PrintMinMax“ of the x-coordinate (longitude) as well as y-coordinate (latitude) look like this:

Variable: shp
Type: file
File path: /home/States_shapefile_data/GADM/TX_OK_KS_CO.nc
Number of global attributes: 11
Number of dimensions: 5
Number of variables: 4
(0) min=0   max=0

ncdump -h TX_OK_KS_CO.nc returns this:

netcdf TX_OK_KS_CO {
dimensions:
geometry = 2 ;
segments = 2 ;
num_features = 4 ;
num_segments = 262 ;
num_points = 102749 ;
variables:
int geometry(num_features, geometry) ;
int segments(num_segments, segments) ;
double x(num_points) ;
double y(num_points) ;

// global attributes:
:segs_numPnts = 1 ;
:segs_xyzIndex = 0 ;
:geom_numSegs = 1 ;
:geom_segIndex = 0 ;
:geometry_type = "polygon" ;
:layer_name = "TX_OK_KS_CO" ;
:creation_date = "Mon May 22 17:20:46 EDT 2017" ;
:NCL_Version = "6.4.0" ;
:system = "Linux ou.edu<http://ou.edu/> x86_64 #1 SMP Fri Mar 3 00:04:05 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux" ;
:Conventions = "None" ;
:title = "NCL: convert-OGR-to-netCDF" ;
}

I am not sure how to use this information, or if this is even the information I should be getting, to define the study region when extracting data from my global forecast dataset. I am fairly sure the shapefile is correct. Opening it again in QGIS showed the four states with their outlines and filling, and nothing else.

I appreciate any help!

Thank you very much,
Toni


Toni Klemm
Ph.D. Candidate
The University of Oklahoma
South Central Climate Science Center
phone: 405 325 1272<tel:(405)%20325-1272>
www.toni-klemm.de<https://urldefense.proofpoint.com/v2/url?u=http-3A__www.toni-2Dklemm.de&d=DwMFaQ&c=qKdtBuuu6dQK9MsRUVJ2DPXW6oayO8fu4TfEHS8sGNk&r=trcSnQNnUT2ddBXcpd358A&m=_2unS1RvDS7egGpkWknJW0V9KMHhv2GWDckbCG7ZOcQ&s=QAW2_xJ9gp71I8wBvMBw6JI7Kb7j7d8KF0tlyNuhEYE&e=>
Early Career Climate Forum<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.eccforum.org&d=DwMFaQ&c=qKdtBuuu6dQK9MsRUVJ2DPXW6oayO8fu4TfEHS8sGNk&r=trcSnQNnUT2ddBXcpd358A&m=_2unS1RvDS7egGpkWknJW0V9KMHhv2GWDckbCG7ZOcQ&s=BsFxGtcxiOomgvIRjvEfgPA7YppuvrO0n-TW0SLSiok&e=>
Twitter: @toniklemm<https://urldefense.proofpoint.com/v2/url?u=https-3A__twitter.com_ToniKlemm&d=DwMFaQ&c=qKdtBuuu6dQK9MsRUVJ2DPXW6oayO8fu4TfEHS8sGNk&r=trcSnQNnUT2ddBXcpd358A&m=_2unS1RvDS7egGpkWknJW0V9KMHhv2GWDckbCG7ZOcQ&s=B5YLwKe0WRmXO4D2nzE_HgiPLcFnIFvPX6Z-kf3KSrw&e=>







_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk<https://urldefense.proofpoint.com/v2/url?u=http-3A__mailman.ucar.edu_mailman_listinfo_ncl-2Dtalk&d=DwMFaQ&c=qKdtBuuu6dQK9MsRUVJ2DPXW6oayO8fu4TfEHS8sGNk&r=trcSnQNnUT2ddBXcpd358A&m=_2unS1RvDS7egGpkWknJW0V9KMHhv2GWDckbCG7ZOcQ&s=g_J7lY2OPhdO6e1uZkyCAbj2G1XmDbXoRk01EanWTdI&e=>



-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170524/f83e492c/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: GFDL_TAS_bias_1980-2014_target_month_01_shp.png
Type: image/png
Size: 103416 bytes
Desc: GFDL_TAS_bias_1980-2014_target_month_01_shp.png
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170524/f83e492c/attachment-0001.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: GFDL_TAS_bias_1980-2014_target_month_01.nc
Type: application/octet-stream
Size: 6116 bytes
Desc: GFDL_TAS_bias_1980-2014_target_month_01.nc
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170524/f83e492c/attachment-0003.obj 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: GFDL_TAS_bias_1980-2014_target_month_01_shp.nc
Type: application/octet-stream
Size: 6116 bytes
Desc: GFDL_TAS_bias_1980-2014_target_month_01_shp.nc
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170524/f83e492c/attachment-0004.obj 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 4_states_merged.nc
Type: application/octet-stream
Size: 13484 bytes
Desc: 4_states_merged.nc
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170524/f83e492c/attachment-0005.obj 


More information about the ncl-talk mailing list