[R-sig-Geo] FW: proj4string error on data.frame generated from NCEP data
Loïc Dutrieux
loic.dutrieux at wur.nl
Fri Apr 8 12:27:51 CEST 2016
>>
>> -----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
More information about the R-sig-Geo
mailing list