r read NetCDF and export as shapefile

Question

I need to read a NetCDF file with R and export each time step as a smoothed polygon shapefile. I have two problems: smoothing the raster and exporting to shapefile with proper projection from the NC file. The output is a regular grid and is not projected.

Here is a sample code:

>NCFileName = MyncFile.nc
NCFile = open.ncdf(NCFileName) 
NCFile 
[1] "file CF_OUTPUT.nc has 6 dimensions:"
[1] "time   Size: 61"
[1] "height   Size: 8"
[1] "lat   Size: 185"
[1] "lon   Size: 64"
[1] "Time   Size: 61"
[1] "DateStrLen   Size: 19"
[1] "------------------------"
[1] "file CF_OUTPUT.nc has 20 variables:"
[1] "float temp[lon,lat,height,time]  Longname:Temperature Missval:1e+30"
[1] "float relh[lon,lat,height,time]  Longname:Relative Humidity Missval:1e+30"
[1] "float airm[lon,lat,height,time]  Longname:Air density Missval:1e+30"
[1] "float z[lon,lat,height,time]  Longname:Layer top altitude Missval:1e+30"
[1] "float ZH[lon,lat,height,time]  Longname:Layer top altitude Missval:1e+30"
[1] "float hlay[lon,lat,height,time]  Longname:Layer top altitude Missval:1e+30"
[1] "float PM10ant[lon,lat,height,time]  Longname:PM10ant Concentration Missval:1e+30"
[1] "float PM10bio[lon,lat,height,time]  Longname:PM10bio Concentration Missval:1e+30"
[1] "float PM10[lon,lat,height,time]  Longname:PM10 Concentration Missval:1e+30"
[1] "float PM25ant[lon,lat,height,time]  Longname:PM25ant Concentration Missval:1e+30"
[1] "float PM25bio[lon,lat,height,time]  Longname:PM25bio Concentration Missval:1e+30"
[1] "float PM25[lon,lat,height,time]  Longname:PM25 Concentration Missval:1e+30"
[1] "float C2H4[lon,lat,height,time]  Longname:C2H4 Concentration Missval:1e+30"
[1] "float CO[lon,lat,height,time]  Longname:CO Concentration Missval:1e+30"
[1] "float SO2[lon,lat,height,time]  Longname:SO2 Concentration Missval:1e+30"
[1] "float NO[lon,lat,height,time]  Longname:NO Concentration Missval:1e+30"
[1] "float NO2[lon,lat,height,time]  Longname:NO2 Concentration Missval:1e+30"
[1] "float O3[lon,lat,height,time]  Longname:O3 Concentration Missval:1e+30"
[1] "char Times[DateStrLen,Time]  Longname:Times Missval:NA"
[1] "float HGT[lon,lat,time]  Longname:Topography Missval:1e+30"

nc.a=get.var.ncdf(NCFile , varid = 'NO2', start=c(1,1,1,1), count=c(-1,-1,1,1))
Pol <- rasterToPolygons(raster(nc.a),dissolve = TRUE)
Pol
class       : SpatialPolygonsDataFrame 
features    : 11829 
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 1
names       :             layer 
min values  : 0.219758316874504 
max values  :  0.84041428565979 
writeOGR(Pol, dsn = getwd(), layer = 'testPol', driver = 'ESRI Shapefile', overwrite_layer = TRUE)

What I get, however, are grided polygons that are not projected.

UPDATE: Following @kakk11 and @RobertH answers, I was able to solve part of the problem. I still get a grid-like polygons, not smoothed. Here is what I did so far: I couldn't extract the variable directly to raster as @RobertH suggested. so I used the 'get.var.ncdf' and then 'raster':

NCFileName = 'MyncFile.nc'
NCFile = open.ncdf(NCFileName)
nc.a = get.var.ncdf(NCFile, varid = 'NO2', start=c(1,1,1,13), count=c(-1,-1,1,1))
nc.a = raster(nc.a)
# put in correct extent:
lat  = NCFile$dim$lat$vals
lon  = NCFile$dim$lon$vals
ExtentLat = range(lat)
ExtentLon = range(lon)
rm(lat,lon)

nc.a = flip(t(nc.a), direction='y')

# Give it lat/lon coords 
extent(nc.a) = c(ExtentLon,ExtentLat)

Then the 'cut' command returns vector, so i used 'ratser:reclassify':

cuts = c(0,5,15,30,50)
classes <- cbind(cuts[1:length(cuts)-1],cuts[2:length(cuts)],cuts[2:length(cuts)])
nc.class <- reclassify(nc.a, classes)

I then used the 'rasterToPolygons' with 'dissolve=TRUE' to create the polygons:

pol <- rasterToPolygons(nc.class, dissolve = TRUE)
# set UTM projection:
WGS84_Projection = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(pol) <- CRS(WGS84_Projection)
writeOGR(pol, dsn = getwd(), layer = 'file' , driver = 'ESRI Shapefile', overwrite_layer = TRUE)

Still, all this creates polygon shapefile with the polygons that are not smooth, which is the main challenge.

Could use some help with this.

Ilik

Answer 1

You first need correctly create a RasterLayer, like this:

 r <- raster('MyncFile.nc', var='NO2')
 # or, to get all time steps at once
 # brick('MyncFile.nc', var='NO2')

You could then generalize (classify) the values using reclassify or cut. For example

 cuts <- seq(0.2, 0.9, 0.1)
 rc <- cut(r, cuts)

Make polygons and save to shapefile

 pol <- rasterToPolygons(rc, dissolve = TRUE)
 shapefile(pol, 'file.shp')