[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