Мне нужно прочитать файл NetCDF с помощью R и экспортировать каждый временной шаг как сглаженный шейп-файл многоугольника. У меня две проблемы: сглаживание растра и экспорт в шейп-файл с правильной проекцией из файла ЧПУ. Результатом является регулярная сетка, которая не прогнозируется.
Вот пример кода:
>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)
Однако я получаю многоугольники с сеткой, которые не проецируются.
ОБНОВЛЕНИЕ: после ответов @ kakk11 и @RobertH я смог частично решить проблему. У меня по-прежнему получаются не сглаженные полигоны в виде сетки. Вот что я сделал до сих пор: я не мог извлечь переменную напрямую в растр, как предложил @RobertH. поэтому я использовал get.var.ncdf, а затем 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)
Затем команда 'cut' возвращает вектор, поэтому я использовал '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)
Затем я использовал rasterToPolygons с disolve = TRUE для создания полигонов:
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)
Тем не менее, все это создает шейп-файл многоугольника с негладкими многоугольниками, что является основной проблемой.
Может понадобиться помощь с этим.
Илик
