gabriela bucini
2013-Sep-28 00:09 UTC
[R] netCDF to raster and projection to UTM: output NA
Dear R colleagues,
I'm facing a problem with the projection of a netCDF file.
My original netCDF file is in lat/lon coordinates and I want to
project it in UTM.
I use the raster package with the function"raster" to open the file
and then "projectRaster" to change projection (I want to maintain a
similar cell resolution as the input netCDF). I obtain a raster of two
cells with NA values. Even if I change the resolution for the
projected output, I always obtain a raster of NA values.
Can anybody help me figuring out where I'm making a mistake?
Thank you very much for your help!
Gabriela
This is my code:
nc2r <-
raster("./../Data/BCCA_daily_original/CMIP5/access1-0.1.rcp45/bcca5/Extraction_pr.nc",
varname = "pr", level=1)
# netCDF info
print(nc2r)
# raster info
nc2r
summary(nc2r)
#define the new projection (UTM)
newproj <- "+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m
+no_defs"
# change the projection maintaining a similar resolution to input image
pr2 <- projectRaster(nc2r, crs=newproj, res=90000)
pr2
summary(pr2)
OUTPUTS:
> print(nc2r)
[1] "File
C:\\Users\\GB\\Dropbox\\Downscaling\\Data\\BCCA_daily_original\\CMIP5\\access1-0.1.rcp45\\bcca5\\Extraction_pr.nc
(NC_FORMAT_CLASSIC):"
[1] ""
[1] " 1 variables (excluding dimension variables):"
[1] " float pr[longitude,latitude,time,projection] "
[1] " typeConversion_op_ncl: double converted to float"
[1] " _FillValue: 1.00000002004088e+20"
[1] " associated_files: baseURL:
http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile:
gridspec_atmos_fx_ACCESS1-0_historical_r0i0p0.nc areacella:
areacella_fx_ACCESS1-0_historical_r0i0p0.nc"
[1] " missing_value: 1.00000002004088e+20"
[1] " history: 2012-01-17T12:34:12Z altered by CMOR:
replaced missing value flag (-1.07374e+09) with standard missing value
(1e+20)."
[1] " cell_measures: area: areacella"
[1] " cell_methods: time: mean"
[1] " units: mm/d"
[1] " comment: at surface; includes both liquid and solid
phases from all types of clouds (both large-scale and convective)"
[1] " long_name: Precipitation"
[1] " standard_name: precipitation_flux"
[1] " time: 0.5"
[1] ""
[1] " 4 dimensions:"
[1] " latitude Size:13"
[1] " valid_max: 52.875"
[1] " long_name: Latitude"
[1] " valid_min: 25.125"
[1] " units: degrees_north"
[1] " axis: Y"
[1] " longitude Size:10"
[1] " long_name: Longitude"
[1] " axis: X"
[1] " units: degrees_east"
[1] " modulo: 360"
[1] " topology: circular"
[1] " time Size:54787"
[1] " calendar: standard"
[1] " units: days since 1950-01-01 00:00:00"
[1] " standard_name: time"
[1] " long_name: time"
[1] " axis: T"
[1] " projection Size:1 *** is unlimited ***"
[1] ""
[1] " 15 global attributes:"
[1] " institution: CSIRO (Commonwealth Scientific and
Industrial Research Organisation, Australia), and BOM (Bureau of
Meteorology, Australia)"
[1] " institute_id: CSIRO-BOM"
[1] " model_id: ACCESS1-0"
[1] " frequency: day"
[1] " experiment: historical"
[1] " experiment_id: historical"
[1] " parent_experiment: pre-industrial control"
[1] " parent_experiment_id: piControl"
[1] " parent_experiment_rip: r1i1p1"
[1] " creation_date: Fri Sep 14 14:51:45 PDT 2012"
[1] " references: Daily BC method: modified version of Maurer
EP, Hidalgo HG, Das T, Dettinger MD, Cayan DR, 2010, Hydrol Earth Syst
Sci 14:1125-1138\nCA method: Hidalgo HG, Dettinger MD, Cayan DR, 2008,
California Energy Commission technical report
CEC-500-2007-123\nReference period obs: updated version of Maurer EP,
Wood AW, Adam JC, Lettenmaier DP, Nijssen B, 2002, J Climate
15(22):3237???3251, \nprovided via
http://www.engr.scu.edu/~emaurer/gridded_obs/index_gridded_obs.html"
[1] " contacts: Bridget Thrasher:
bridget at climateanalyticsgroup.org or Ed Maurer: emaurer at scu.edu"
[1] " documentation: http://gdo-dcp.ucllnl.org"
[1] " NCO: 4.0.8"
[1] " Projections: access1-0.1.rcp45, "
# raster info
> nc2r
class : RasterLayer
dimensions : 13, 10, 130 (nrow, ncol, ncell)
resolution : 0.125, 0.125 (x, y)
extent : 286.625, 287.875, 43.875, 45.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Gabriela
Bucini\Dropbox\Downscaling\Data\BCCA_daily_original\CMIP5\access1-0.1.rcp45\bcca5\Extraction_pr.nc
names : Precipitation
z-value : 1
zvar : pr
level : 1
class : RasterLayer
dimensions : 13, 10, 130 (nrow, ncol, ncell)
resolution : 0.125, 0.125 (x, y)
extent : 286.625, 287.875, 43.875, 45.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Gabriela
Bucini\Dropbox\Downscaling\Data\BCCA_daily_original\CMIP5\access1-0.1.rcp45\bcca5\Extraction_pr.nc
names : Precipitation
z-value : 1
zvar : pr
level : 1
> summary(values(nc2r))
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.5352 1.2530 1.6820 1.7590 2.2790 3.4180 6
#define the new projection
newproj <- "+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m
+no_defs"
# change the projection maintaining a similar resolution to input image
pr2 <- projectRaster(nc2r, crs=newproj, res= 90000)
> pr2
class : RasterLayer
dimensions : 2, 1, 2 (nrow, ncol, ncell)
resolution : 90000, 90000 (x, y)
extent : 626955.8, 716955.8, 4862518, 5042518 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m
+no_defs +towgs84=0,0,0
data source : in memory
names : Precipitation
values : NA, NA (min, max)
> summary(pr2)
Precipitation
Min. NA
1st Qu. NA
Median NA
3rd Qu. NA
Max. NA
NA's 2