This repository has been archived by the owner on Nov 27, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 1
/
FromNCDFSG.R
122 lines (103 loc) · 3.78 KB
/
FromNCDFSG.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#'@title Convert NetCDF to sp objects
#'
#'
#'@param nc_file A string file path to the nc file to be read.
#'
#'@description
#'Attemps to convert a NetCDF-CF DSG Simple Geometry file into an sp object.
#'
#'@references
#'https://github.com/bekozi/netCDF-CF-simple-geometry
#'
#'@importFrom ncdf4 nc_open nc_close ncvar_get ncatt_get
#'@importFrom sp Polygon Polygons SpatialPolygons SpatialPolygonsDataFrame CRS Line Lines SpatialLines SpatialLinesDataFrame SpatialPointsDataFrame
#'@importFrom netcdf.dsg read_instance_data
#'
#'@export
FromNCDFSG = function(nc_file) {
nc <- nc_open(nc_file)
checkVals <- checkNCDF(nc)
instance_id<-checkVals$instance_id
instanceDim<-checkVals$instanceDim
geom_container <- checkVals$geom_container
variable_list <- checkVals$variable_list
crs <- checkVals$crs
line<-FALSE; poly<-FALSE; point<-FALSE
if(grepl("polygon", geom_container$geom_type)) { poly<-TRUE
} else if(grepl("line", geom_container$geom_type)) { line<-TRUE
} else point <- TRUE
xCoords <- c(ncvar_get(nc, geom_container$x))
yCoords <- c(ncvar_get(nc, geom_container$y))
if(length(crs) == 0) {
prj <- "+proj=longlat +datum=WGS84"
} else {
prj <- getPrjFromNCDF(crs)
}
if(point) {
point_data <- matrix(c(xCoords,
yCoords), ncol=2)
dataFrame <- read_instance_data(nc, instanceDim)
if(geom_container$geom_type == "multipoint") {
stop("reading multipoint is not supported yet.")
# This is where handling for multipoint would go.
}
SPGeom <- SpatialPointsDataFrame(point_data, proj4string = CRS(prj),
data = dataFrame, match.ID = FALSE)
} else {
node_count <- c(ncvar_get(nc, geom_container$node_count))
if(!is.null(instance_id)) {
instance_names <- ncvar_get(nc, instance_id)
} else {
instance_names <- as.character(c(1:length(node_count)))
}
if(is.character(geom_container$part_node_count)) {
part_node_count <- ncvar_get(nc, geom_container$part_node_count)
} else {
part_node_count <- node_count
}
if(is.character(geom_container$part_type)) {
part_type <- ncvar_get(nc, geom_container$part_type)
} else {
part_type <- rep(pkg.env$multi_val, length(part_node_count))
}
node_start <- 1
geom_node_stop <- 0
pInd <- 1
Srl <- list()
for(geom in 1:length(node_count)) {
geom_node_stop <- geom_node_stop + node_count[geom]
srl <- list()
while(node_start < geom_node_stop) {
part_node_stop <- node_start + part_node_count[pInd] - 1
if(part_type[pInd] == pkg.env$hole_val) { hole <- TRUE
} else { hole <- FALSE }
coords <- matrix(c(xCoords[node_start:part_node_stop],yCoords[node_start:part_node_stop]),ncol=2)
if(poly) { tsrl<-Polygon(coords, hole=hole)
} else if(line) { tsrl<-Line(coords) }
dimnames(tsrl@coords) <- list(NULL, c("x", "y"))
srl <- append(srl, tsrl)
node_start <- node_start + part_node_count[pInd]; pInd <- pInd + 1
}
if(poly) {
Srl <- append(Srl, Polygons(srl, as.character(geom)))
} else if(line) {
Srl <- append(Srl, Lines(srl, as.character(geom)))
}
}
dataFrame <- read_instance_data(nc, instanceDim)
for(varName in names(dataFrame)) {
if(!varName %in% variable_list) {
dataFrame[varName] <- NULL
}
}
if(poly) {
SPGeom <- SpatialPolygonsDataFrame(SpatialPolygons(Srl, proj4string = CRS("+proj=longlat +datum=WGS84")),
dataFrame, match.ID = FALSE)
} else if(line) {
SPGeom <- SpatialLinesDataFrame(SpatialLines(Srl, proj4string = CRS("+proj=longlat +datum=WGS84")),
dataFrame, match.ID = FALSE)
}
}
nc_close(nc)
return(SPGeom)
}