1

NARCCAPの温度 .nc ファイルを使用しています。これらのデータには極平射図法があります。気温の最低値と最高値から、メープル シロップの生産日と見なされる日数のマトリックスを作成しました。

このマトリックスをラスターに変換し、このラスターを経度/緯度投影に投影したいと考えています。

## This is the metadata for the projection from the .nc file:

 # float lat[xc,yc]   
 #            long_name: latitude
 #            standard_name: latitude
 #            units: degrees_north
 #            axis: Y
 #  float lon[xc,yc]   
 #            long_name: longitude
 #            standard_name: longitude
 #            units: degrees_east
 #            axis: X
# float tasmax[xc,yc,time]   
#             coordinates: lon lat level
#             _FillValue: 1.00000002004088e+20
#             original_units: K
#             long_name: Maximum Daily Surface Air Temperature
#             missing_value: 1.00000002004088e+20
#             original_name: T_MAX_GDS5_HTGL
#             units: K
#             standard_name: air_temperature
#             cell_methods: time: maximum (interval: 24 hours)
#             grid_mapping: polar_stereographic

# grid_mapping_name: polar_stereographic
# latitude_of_projection_origin: 90
# standard_parallel: 60
# false_easting: 4700000
# false_northing: 8400000
# longitude_of_central_meridian: 263
# straight_vertical_longitude_from_pole: 263

# The production days matrix I've created is called from a saved file:
path.ecp2 <- paste0("E:/all_files/production/narccap/GFDL/Production_Days_SkinnerECP2", 
               year, ".RData")
file.ecp2 <- get(load(path.ecp2))
dim(file.ecp2)
# 147 116
rast.ecp2 <- raster(file.ecp2)
rast.ecp2 <- flip(t(rast.ecp2), 2)
# class       : RasterLayer 
# dimensions  : 116, 147, 17052  (nrow, ncol, ncell)
# resolution  : 0.006802721, 0.00862069  (x, y)
# extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
# coord. ref. : NA 
# data source : in memory
# names       : layer 
# values      : 0, 671  (min, max)

# I assign the polar stereographic crs to this production days raster:
crs("+init=epsg:3031")
ecp2.proj <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=4700000 +y_0=8400000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
crs(rast.ecp2) <- crs(ecp2.proj)

rast.ecp2
# class       : RasterLayer 
# dimensions  : 116, 147, 17052  (nrow, ncol, ncell)
# resolution  : 0.006802721, 0.00862069  (x, y)
# extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=4700000 +y_0=8400000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : layer 
# values      : 0, 671  (min, max)

以前に機能した手順を使用すると (こちらを参照)、rast.ecp2 の値はすべて NA になります。どこが間違っていますか？

# The projection I want to project TO:
source_rast <- raster(nrow=222, ncol=462, xmn=-124.75, xmx=-67, ymn=25.125, ymx=52.875,
                      crs="+proj=longlat +datum=WGS84")
rast.ecp2LL <- projectRaster(rast.ecp2, source_rast)

rast.ecp2LL
# class       : RasterLayer 
# dimensions  : 222, 462, 102564  (nrow, ncol, ncell)
# resolution  : 0.125, 0.125  (x, y)
# extent      : -124.75, -67, 25.125, 52.875  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : layer 
# values      : NA, NA  (min, max)
4

1 に答える 1

1

私が見つけた解決策を投稿しています。これは、この投稿と回答に基づいています。まず、.nc ファイルの xc 座標と yc 座標を経度と緯度のポイントに変換する必要がありました。その後、ラスターを適切に再投影できました。以下は動作したコードです。

mycrs.nc ファイルに「付属している」CRS であることに注意してください。SpatialPointsxc/yc から変換しSpatialPointsて関連する CRS をドロップするため、 に割り当てる必要があります。

years <- seq(from=1971, to=2000, by=5)
model <- "CRCM"

convert.lonlat <- function(model, year) 
{
  max.stem <- "E:/all_files/www.earthsystemgrid.org/CCSM/tasmax_"
  inputfile <- paste0(max.stem, model, "_ccsm_", year, "010106.nc")
  lat <- raster(inputfile, varname="lat")
  lon <- raster(inputfile, varname = "lon")
  plat <- rasterToPoints(lat)
  plon <- rasterToPoints(lon)
  lonlat <- cbind(plon[,3], plat[,3])
  lonlat <- SpatialPoints(lonlat, proj4string = crs(base.proj))
  mycrs <- crs("+proj=stere +lon_0=263 +x_0=3475000 +y_0=7475000 +lat_0=90 +ellps=WGS84")
  plonlat <- spTransform(lonlat, CRSobj = mycrs)
  maxs <- brick(inputfile, varname="tasmax")
  projection(maxs) <- mycrs
  extent(maxs) <- extent(plonlat)
  max.lonlat <- projectRaster(maxs, base.proj)
  save(max.lonlat, file=paste0("E:/all_files/production/narccap/CCSM/", model, "max_lonlat_", year, ".RData"))

  min.stem <- "E:/all_files/www.earthsystemgrid.org/CCSM/tasmin_"
  inputfile <- paste0(min.stem, model, "_ccsm_", year, "010106.nc")
  lat <- raster(inputfile, varname="lat")
  lon <- raster(inputfile, varname = "lon")
  plat <- rasterToPoints(lat)
  plon <- rasterToPoints(lon)
  lonlat <- cbind(plon[,3], plat[,3])
  lonlat <- SpatialPoints(lonlat, proj4string = crs(maurer.proj))
  mycrs <- crs("+proj=stere +lon_0=263 +x_0=3475000 +y_0=7475000 +lat_0=90 +ellps=WGS84")
  plonlat <- spTransform(lonlat, CRSobj = mycrs)
  mins <- brick(inputfile, varname="tasmin")
  projection(mins) <- mycrs
  extent(mins) <- extent(plonlat)
  min.lonlat <- projectRaster(mins, maurer.proj)
  save(min.lonlat, file=paste0("E:/all_files/production/narccap/CCSM/", model, "min_lonlat_", year, ".RData"))
}

lapply(years, convert.lonlat, model=model)

ここから、保存したファイルmax.lonlatとmin.lonlat.

ありがとう！これが役に立てば幸いです。

于 2016-08-15T17:51:12.580 に答える