[R-sig-Geo] proj4string error on data.frame generated from NCEP data

Maggie CY Lau maglau at Princeton.EDU
Sat Apr 9 13:34:16 CEST 2016


Frede, Loïc, Edzer, thank you very much for your responses.

The coding lines provided Frede and Loïc worked. The tif file was imported into ArcMap and could be viewed in different projections. 


On Apr 8, 2016, at 6:27 AM, Loïc Dutrieux wrote:

> 
>>> 
>>> -----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
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Maggie C.Y. Lau, PhD
Department of Geosciences
B80 Guyot Hall
Princeton University
Princeton, NJ 08544, US
Cell: +1-609-356-8145



More information about the R-sig-Geo mailing list