I am analyzing seasonal climate forecast skill (netCDF data) for a region comprised of the U.S. states of Texas, Oklahoma, Kansas, and Colorado using NCL (NCAR Command Language). Can I use shapefiles of these states to create a spatial subset of the data so it covers only these states, and if so, how? At the moment I use latitudes and longitudes to define my study region (see the example plot attached), but some of the most extreme model outliers are outside my study region (for example, Mexico, Missouri, Iowa), which distorts my results. Can anyone help?
Thank you very much!
Hi Toni. You can use the 'mask' function of ncl for this purpose. Refer to the attached snapshots that I got from the user guide. Track the usage of the file stored in the variable 'file1' in the script there. However to store your masking file like that, first you need to convert it to .nc file. You can easily do that using ncl_convert 2nc function.
Hope this helps. Cheers! Pallav
I think you should convert your data into preferable format then you can mask it through your shapefile.
Hi Pallav! It sounds like that's exactly what I was looking for. Thanks a lot!
How do I get the mask file? Are there templates for often used areas in common grid cell resolutions (0.5 x 0.5 degrees is what I need), or do I need to create it myself?
Thanks again!
Avinash, thanks! My data is already in netCDF (.nc) format. It seems that's what I need to use the mask function.
Toni. Yes you gotta create it yourself. It isn't hard at all.
You have your data (forecast skill) in nc. I understand that you already have the shape for of your study area (particular States). Thus the mask is in shp. The "ncl_convert2nc" can be used to convert your shp file to nc mask file. BTW ncl_convert2nc is one of the tool/ function of NCL itself. However, strangely I got only one hit for the tool in the ncl user guide (attached pic). You should be able to find it's syntax and usage examples over the Internet though.
All the best!
Yes Pallav, I think so. I have a shapefile with all US states (from NCAR: http://www.ncl.ucar.edu/Applications/Data/#shp --> states.shp) and I managed to convert it from .shp to .nc by running "ncl_convert2nc states.shp states.nc -nc2c" directly from the command line (not writing an NCL script).
However, I'm not sure I understand what I am getting here and how to use it. I attached a screenshot of "ncdump -h states.nc". Many of the variables listed are just lists of numbers, without context or meaning. How do I tell NCL what my study region is?
Toni. the solution is simple. I bet you have a good grip on GIS.
1. Using the .shp file that you have, extract the polygons corresponding to the "states of interest" in GIS. In ArcGIS you can do that using "copy features" tool.
2. Now convert this new .shp file (say states_Toni.shp) to .nc (say States_Toni.nc) using NCL as you did before.
3. Finally mask your forecast skill .nc file using States_Toni.nc. That should be it.
I hope you succeed on this. Cheers!
Thanks again, Pallav! It turns out that you can't use QGIS to create a shapefile that works with ncl_convert2nc. When converting the QGIS shapefile, for some reason all x- and y-coordinates had values of zero which of course makes it unusable.
I had to use ArcGIS to create the shapefile. After creating the shapefile in ArcMap, the conversion to .nc worked fine, and I'm using this file now to run over my data and create a subset based on the shapefile coordinates! There are a few examples on the NCL website. This one worked for me for both creating a data subset and for plotting the subset: https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl
You can easily do this in R studio using raster and netcdf4 packages.
library(raster)
library(ncdf4)
library(maptools)
poly_shape=readShapePoly("yyyyyyy.shp")#shapefile name
data=brick("xxxxxxxxxx.nc") #nc file name
output=mask(data, poly_shape)
* note:make sure that projection of nc and shape files is same
Amrit Thapa
Thanks for the code. I tried using this but my netCDF file has a mercator projection and that doesn't let me mask the shapefile over netCDF appropriately. Do you have a fix for this?
This can be done in R using the raster package.
https://cran.r-project.org/web/packages/raster/index.html
# R Code ---------------------------------------------------------------------------------------------------------------
library(raster)
# Read in shapefile and netCDF
shp = shapefile('C:/path/to/shapefile/ILshapefile.shp')
br = brick('C:/path/to/netcdf/hurs_bced_1960_1999_gfdl-esm2m_hist_1950-2004.nc4')
# Transform shapefile to match the CRS of the netCDF
shp = spTransform(shp, crs(br))
# Mask all cells outside of shapefile
br.subset = mask(br, shp)
# See the result
plot(br.subset, 1)
plot(shp, add=T)
##copied from Jacob Fraser
- Badges
- Science topic
- Similar topics
- Correction