proj4string error on data.frame generated from NCEP data
-----Original Message----- From: Frede Aakmann T?gersen Sent: 8. april 2016 08:01 To: 'Maggie CY Lau'; r-sig-geo at r-project.org Subject: RE: proj4string error on data.frame generated from NCEP data Hi Maggie Your error message says that there is some coordinates with latitudes above 90 degrees. The problem seems to arise from the call to gridded. See below for a dput'ed version of the dataframe "air.sig995.ag.df" called "friend" in my code. This is what I found out:
coordinates(friend) <- ~longitude+latitude bbox(friend)
min max longitude -177.5 180 latitude 60.0 90
apply(coordinates(friend), 2, range)
longitude latitude [1,] -177.5 60 [2,] 180.0 90 After using gridded() the bounding box is enlarged somewhat but the coordinates still unchanged.
gridded(friend) <- TRUE bbox(friend)
min max longitude -178.75 181.25 latitude 58.75 91.25
apply(coordinates(friend), 2, range)
longitude latitude [1,] -177.5 60 [2,] 180.0 90
If you leave out the gridded() part I guess you can proceed towards your goal. Here I did: ## this was okay
crs <- "+proj=longlat +datum=WGS84" proj4string(friend) <- crs
## see attached plot
png("longlat.png")
plot(friend)
dev.off()
Do a reprojection with spTransform():
stere <- "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" friend.stere <- spTransform(friend, CRS(stere))
And the plot is attached:
png("stere.png")
plot(friend.stere)
dev.off()
X11cairo Why gridded() is extending domain outside permissible values for longitude and latitude I don't know.
My guess is that gridded assumes coordinates provided correspond to pixels centers, while top-left coordinates are provided.
However, I don't understand why raster does not complain about a raster
having an extent beyond the boundaries of its CRS.
If I do:
library(RNCEP)
library(raster)
air.sig995.extent <- NCEP.gather(variable='air.sig995', level='surface',
months.minmax=c(7),
years.minmax=c(2011),
lat.southnorth=c(60,89.9), lon.westeast=c(-179,180),
reanalysis2=FALSE, return.units=TRUE)
air.sig995.ag <- NCEP.aggregate(wx.data=air.sig995.extent, YEARS=TRUE,
MONTHS=TRUE,
DAYS=FALSE, HOURS=FALSE, fxn='mean')
air.sig995.ag.df <- NCEP.array2df(air.sig995.ag[,,1], var.names =
'Temperature')
# Reorder columns
df <- air.sig995.ag.df[,c(2,1,3)]
r <- rasterFromXYZ(df, crs = CRS('+proj=longlat +datum=WGS84'))
r
plot(r)
writeRaster(r, '~/sandbox/arctic.tif')
everything works... (but is probably wrong)
Cheers,
Lo?c
Cheers, Lo?c
\Frede ##### friend is a copy of air.sig995.ag.df <- NCEP.array2df(air.sig995.ag[,,1], var.names = 'Temperature') ## DELETED Yours sincerely / Med venlig hilsen Frede Aakmann T?gersen Specialist, M.Sc., Ph.D. Plant Performance & Modeling Technology & Service Solutions T +45 9730 5135 M +45 2547 6050 frtog at vestas.com http://www.vestas.com Company reg. name: Vestas Wind Systems A/S This e-mail is subject to our e-mail disclaimer statement. Please refer to www.vestas.com/legal/notice If you have received this e-mail in error please contact the sender. -----Original Message----- From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Maggie CY Lau Sent: 7. april 2016 16:24 To: r-sig-geo at r-project.org Subject: [R-sig-Geo] proj4string error on data.frame generated from NCEP data HI all, I am an ecologist who would like to use several datasets generated by NCEP reanalysis. I attempted to use the R package "RNCEP" to extract the data and generate .tiff files. These files will be analyzed by a student using ArcMap. I created the data fame with the information I need, but had problem exporting the data into .tiff file that can be projected in polar stereograph in ArcMap. The use of "+proj=longlat + datum=WGS84" gave the following error message: Error in `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) : Geographical CRS given to non-conformant data: 91.25 Then I tried using ""+proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs". proj4string went fine but the export tiff appeared as a box image instead of an angular one. Please kindly check the following codes from the line "library(sp)" onward. Thanks, Maggie Lau ==== install.packages("RNCEP", dependencies=TRUE) install.packages(c('abind', 'maps', 'fields', 'tgp', 'fossil'), dependencies=TRUE) install.packages("raster") install.packages("rgdal") library(RNCEP) air.sig995.extent <- NCEP.gather(variable='air.sig995', level='surface', months.minmax=c(7), years.minmax=c(2011), lat.southnorth=c(60,89.9), lon.westeast=c(-179,180), reanalysis2=FALSE, return.units=TRUE) air.sig995.ag <- NCEP.aggregate(wx.data=air.sig995.extent, YEARS=TRUE, MONTHS=TRUE, DAYS=FALSE, HOURS=FALSE, fxn='mean') dimnames(air.sig995.ag)[[3]][1] air.sig995.ag.df <- NCEP.array2df(air.sig995.ag[,,1], var.names = 'Temperature') library(sp) coordinates(air.sig995.ag.df) <- ~longitude+latitude gridded(air.sig995.ag.df) <- TRUE crs <- "+proj=longlat + datum=WGS84" proj4string(air.sig995.ag.df) <- CRS(crs) #stere <- "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 # +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" #proj4string(air.sig995.ag.df) <- CRS(stere) library(raster) ras <- raster(air.sig995.ag.df) writeRaster(ras, 'air.sig995.Jul.ag.tif', overwrite=TRUE) [[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo