Converting .shp to .nc
Added by Josh Butterworth over 5 years ago
Hi guys,
I've been reading through an old post by Karin showing how you can use NCL to convert a .shp file to .nc (below). This is to create a netcdf file of marine protected areas in the North Sea for which I have a shapefile.
;-- copy shapefile_utils.ncl from
;-- https://www.ncl.ucar.edu/Applications/Scripts/shapefile_utils.ncl
load "$HOME/NCL/shapefiles/shapefile_utils.ncl"
shpname = "./India/india.shp"
maskname = "./India/india_mask.nc"
print_shapefile_info(shpname)
;-- open data file to get the grid to be used
f = addfile("monmean2007_celcius.nc","r")
;-- read variable
var = f->T2(0,:,:)
var@lat2d = f->XLAT
var@lon2d = f->XLONG
;-- shapefile mask resources
opt = True
opt@return_mask = True ;-- this forces the return of a 0s and 1s mask array
;-- create the mask based on the given shapefile
mask_array = shapefile_mask_data(var, shpname, opt)
mask_array!0 = "y"
mask_array!1 = "x"
mask_array@coordinates = "XLAT XLONG"
;-- create new netCDF file and write mask array
system("rm -f " + maskname)
fout = addfile(maskname,"c")
fout->mask_array = mask_array
fout->XLAT = f->XLAT
fout->XLONG = f->XLONG
Using my data (attached) I've managed to get to the line below by substituting the T2 with time, XLAT with lat and XLONG with lon.
mask_array = shapefile_mask_data(var,shpname,opt)
But then I get the following error message:
(0) shapefile_mask_data: Error: not a valid rectilinear, curvilinear, or uns
tructured grid
Does anyone know why this is happening and where I've gone wrong?
Many thanks,
Josh
SST.nc (204 KB) SST.nc | Data with desired grid | ||
shapefiles.tar (7.83 MB) shapefiles.tar | Shapefiles |
Replies (3)
RE: Converting .shp to .nc - Added by Josh Butterworth over 5 years ago
I realise in the file I uploaded after removing the time axis to enable the upload I no longer have a 'time' var, can this be done without including 'time'? If not I've attached the same data but with a years worth of time values.
RE: Converting .shp to .nc - Added by Karin Meier-Fleischer over 5 years ago
Hi Josh,
this is an NCL question but I've seen that you have send it to ncl-talk too. S, you will get two answers.
It doesn't matter if you have a time dimension. The adaption of the NCL script above gives a netCDF file of the masked variable.
load "$HOME/NCL/shapefiles/shapefile_utils.ncl" shpname = "ospar_polygon_wdpa_simplified.shp" maskname = "ospar_polygon_wdpa_simplified.nc" print_shapefile_info(shpname) ;-- open data file to get the grid to be used f = addfile("SST.nc","r") ;-- read variable var = f->SST(0,:,:) var@lat1d = f->LAT var@lon1d = f->LON ;-- shapefile mask resources opt = True opt@return_mask = True ;-- this forces the return of a 0s and 1s mask array ;-- create the mask based on the given shapefile mask_array = shapefile_mask_data(var, shpname, opt) mask_array!0 = "y" mask_array!1 = "x" mask_array@coordinates = "LAT LON" ;-- create new netCDF file and write mask array system("rm -f " + maskname) fout = addfile(maskname,"c") fout->mask_array = mask_array fout->LAT = f->LAT fout->LON = f->LON
ncdump -h ospar_polygon_wdpa_simplified.nc
netcdf ospar_polygon_wdpa_simplified { dimensions: y = 121 ; x = 211 ; variables: int mask_array(y, x) ; mask_array:coordinates = "LAT LON" ; mask_array:_FillValue = -2147483647 ; } (ncl6.6) [~/Downloads] > ncdump -h ospar_polygon_wdpa_simplified.nc netcdf ospar_polygon_wdpa_simplified { dimensions: y = 121 ; x = 211 ; LAT = 121 ; LON = 211 ; variables: int mask_array(y, x) ; mask_array:coordinates = "LAT LON" ; mask_array:_FillValue = -2147483647 ; float LAT(LAT) ; LAT:long_name = "Latitude" ; LAT:units = "degrees_north" ; LAT:axis = "Y" ; LAT:bounds = "LAT_bnds" ; LAT:point_spacing = "uneven" ; LAT:standard_name = "latitude" ; float LON(LON) ; LON:long_name = "Longitude" ; LON:units = "degrees_east" ; LON:axis = "X" ; LON:bounds = "LON_bnds" ; LON:modulo = 360.f ; LON:point_spacing = "uneven" ; LON:standard_name = "longitude" ; }
-Karin
RE: Converting .shp to .nc - Added by Dana kh over 2 years ago
why the mask.nc after converting the type of lat and lon 2D not Geo2d
how to fix it ???