Dear Robert,

Thanks very much. I adjusted my code to what you sent me and I still have
the problem so I now know that it is with the data although I dont exactly
know what. Perhaps I do not have sufficient variation. I am thinking of
using the altitude data I have as well. Do you think there is any way of
incorporating altitude into this?

Thanks

minti

On Thu, Oct 11, 2012 at 6:02 AM, Robert J. Hijmans <r.hijmans@gmail.com>wrote:

>
> Minti,
> You approach seems fine in general. Below is what I might do, but it is
> essentially the same (there is no role for rasterize). The numbers seem to
> make sense in this case. I cannot comment on why it does not work with your
> real data, but perhaps you can work your way back from here with this
> example.
> Robert
>
>
> x <- 1:10
> y <- 1:10
> z <- 1:10
> d <- data.frame(x,y,z)
>
> library(raster)
> xy <- cbind(x, y)
> # get a reasonable extent
> e <- extent( as.vector(apply(xy, 2, range)) * 1.25 )
> # and a raster with that extent and an arbitrary resolution
> r <- raster(ncols=50, nrows=50, ext=e)
>
> library(gstat)
> IDW <- gstat(id = "z", formula = z~1, locations = ~x+y, data=d,
> set=list(idp =2))
> IDWR <- interpolate(r, IDW)
>
> plot(IDWR)
> points(d[, 1:2], cex=d[,3]/2)
>
>
>
> On Wed, Oct 10, 2012 at 1:34 PM, alemante zd <mintibezabih@gmail.com>wrote:
>
>> Dear Robert,
>>
>> What I am trying to do is interpolate a station level temprature data to
>> farm level and I basically used the code below. The problem was that my
>> interpolated results seemed very implausible (for instance if my
>> temperature observations ranged between 6 and 28, I ended up getting
>> interpolated result between 15 and 17 and extreme temperatures are critical
>> to statistical analysis I will be doing after the interpolation). So I
>> thought the reason I was getting these results was because somehow my
>> raster formation was bad. . Perhaps I have a different problem than the
>> raster, but I have tried different interpolation mathods (thin plate
>> spline, kriging and IDW) but my problem will not go away. Anyway here is my
>> full code:
>>
>> So what I try to do is create a raster from my farm level location
>> observations and then interpolate using station level location and
>> temperature observations
>>
>> cheers
>> minti
>>
>> # In my actual code, I would read a file containing the latititude and
>> longitude of my farm locations
>> x <- 1:10
>> y <- 1:10
>>
>>
>>
>> F1<-data.frame(x, y)
>> library(raster)
>> r <- raster(ncols=36, nrows=18)
>> xy <- cbind(x, y)
>> r1 <- rasterize(xy, r)
>>
>> # In my actual code, I would read a file containing the temperature
>> observation and latititude and longitude of my stations
>>
>> x <- 1:10
>>  y <- 1:10
>> z<- 1:10
>>  da1<-data.frame(x,y,z)
>>
>>
>> IDW <- gstat(id = "z", formula = z~1, locations = ~x+y, data=da1,
>> set=list(idp =2))
>> IDWRaster <- interpolate(r1, IDW)
>>
>>
>> On Wed, Oct 10, 2012 at 5:12 PM, Robert J. Hijmans <r.hijmans@gmail.com>wrote:
>>
>>> Minti,
>>> I do not know what you are trying to achieve so I cannot tell you how to
>>> fix this.
>>> I think we need a clear description of your objectives before it is
>>> possible to give useful advice.
>>> Robert
>>>
>>> On Mon, Oct 8, 2012 at 3:21 PM, alemante zd <mintibezabih@gmail.com>wrote:
>>>
>>>> Dear Robert,
>>>>
>>>> Thanks.  I have now tried rasterizing as per your suggestion but I was
>>>> not sure what should exactly get into the rasterize as the raster.object
>>>> argument (r in my case). Below is a reproducable example I did with what I
>>>> think could be potential r's.
>>>>
>>>>  x <- 1:10
>>>> y <- 1:10
>>>> F1<-data.frame(x, y)
>>>> library(raster)
>>>> r <- raster(as.matrix(F1))
>>>> #r <- raster(ncols=36, nrows=18)
>>>> xy <- cbind(x, y)
>>>> r1 <- rasterize(xy, r)
>>>>
>>>> I would be happy if I get suggestions on how I can fix this.
>>>>
>>>> Thanks
>>>> minti
>>>>  On Mon, Oct 8, 2012 at 4:13 AM, Robert J. Hijmans <r.hijmans@gmail.com
>>>> > wrote:
>>>>
>>>>>
>>>>>
>>>>>  On Sat, Oct 6, 2012 at 3:08 PM, alemante zd <mintibezabih@gmail.com>wrote:
>>>>>
>>>>>> Dear Robert, Tom and Caleb,
>>>>>>
>>>>>> Thanks for your suggestions. I have now fixed my CRs. But I now face
>>>>>> further problems further down in my code as my ultimate objective is
>>>>>> to get a raster layer.
>>>>>>
>>>>>> Here is the full code based on Robert's reproducable example:
>>>>>>
>>>>>>
>>>>>>
>>>>>> x <- 1:10
>>>>>>
>>>>>> y <- 1:10
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> F1<-data.frame(x, y)
>>>>>>
>>>>>> library(sp)
>>>>>>
>>>>>>
>>>>>> #SP <- SpatialPoints(F1,proj4string=CRS("+proj=longlat"))
>>>>>>
>>>>>> SP <- SpatialPoints(F1,proj4string=CRS("+proj=latlong
>>>>>> +datum=potsdam")),
>>>>>>          CRS("+proj=merc"))
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> library(rgdal)
>>>>>>
>>>>>> #spt <- spTransform(SP, CRS("  fill in  the proj4 string you want  ")
>>>>>> )
>>>>>> spt<-spTransform(SpatialPoints(data.frame(x=1:10, y=1:10),
>>>>>>          proj4string=CRS("+proj=latlong +datum=potsdam")),
>>>>>>          CRS("+proj=merc"))
>>>>>>
>>>>>>
>>>>>> gridded(spt) = TRUE
>>>>>>
>>>>>> library(raster)
>>>>>> r = raster(spt)
>>>>>> projection(r) = CRS("+proj=latlong +datum=potsdam")
>>>>>> plot(r)
>>>>>>
>>>>>>
>>>>>> I now get error messages like coordinate intervals not constant and no
>>>>>> vallues associated with the raster layer.
>>>>>>
>>>>>>
>>>>>
>>>>> If you report a problem. Please provide the exact messages you get,
>>>>> not something "like".
>>>>> Either way, the situation seems clear. You have points that are not on
>>>>> a regular grid. To create a raster you could use the raster function
>>>>> rasterize.
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>> Thanks
>>>>>> minti
>>>>>>
>>>>>> On Sat, Oct 6, 2012 at 6:19 PM, Robert J. Hijmans <
>>>>>> r.hijmans@gmail.com> wrote:
>>>>>> >
>>>>>> > Alemante,
>>>>>> >
>>>>>> >> but I keep on getting the message that CRS is wrong.
>>>>>> >> Coordinate Reference System (CRS) arguments: NA
>>>>>> >
>>>>>> > Please show more complete code. Which object are you talking about?
>>>>>> Some of
>>>>>> > the object you create have a CRS (SP) others do not.(pts, F1),
>>>>>> which is
>>>>>> > pretty obvious from reading the sequence of events in your code.
>>>>>> You set the
>>>>>> > CRS for SP, but then you use F1 again which does not have a CRS and
>>>>>> > therefore neither do the objects derived from it. You could do
>>>>>> after,
>>>>>> > gridded(F1) = TRUE,  proj4string(F1) <- CRS("+proj=longlat")
>>>>>> >
>>>>>> >
>>>>>> >> Could you please advice me on how I can get it right?
>>>>>> >
>>>>>> > I need to guess, but I think you want to transform points, in that
>>>>>> case, do
>>>>>> > this:
>>>>>> >
>>>>>> > # example data
>>>>>> > x <- 1:10
>>>>>> > y <- 1:10
>>>>>> >
>>>>>> > F1<-data.frame(x, y)
>>>>>> > library(sp)
>>>>>> > SP <- SpatialPoints(F1,proj4string=CRS("+proj=longlat"))
>>>>>> >
>>>>>> > library(rgdal)
>>>>>> > spt <- spTransform(SP, CRS("  fill in  the proj4 string you want
>>>>>>  ") )
>>>>>> >
>>>>>> >
>>>>>> > Best, Robert
>>>>>> >
>>>>>> >
>>>>>> > On Sat, Oct 6, 2012 at 5:48 AM, alemante zd <mintibezabih@gmail.com>
>>>>>> wrote:
>>>>>> >>
>>>>>> >> Dear all,
>>>>>> >>
>>>>>> >> I would like to sp transform my latitude/longitude vector to by
>>>>>> first
>>>>>> >> telling it what CRS it is. This is what I did.
>>>>>> >>
>>>>>> >> F1<-data.frame(x, y)
>>>>>> >> coordinates(F1)=~x+y
>>>>>> >> SP<-SpatialPoints(F1,proj4string=CRS("+proj=longlat"))
>>>>>> >> pts = SpatialPixels(SpatialPoints(F1))
>>>>>> >> gridded(F1) = TRUE
>>>>>> >>
>>>>>> >> but I keep on getting the message that CRS is wrong.
>>>>>> >> Coordinate Reference System (CRS) arguments: NA
>>>>>> >>
>>>>>> >> Could you please advice me on how I can get it right?
>>>>>> >>
>>>>>> >> cheers
>>>>>> >> minti
>>>>>> >>
>>>>>> >> _______________________________________________
>>>>>> >> R-sig-Geo mailing list
>>>>>> >> R-sig-Geo@r-project.org
>>>>>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>> >
>>>>>> >
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>
>

	[[alternative HTML version deleted]]

