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]]