[R-sig-Geo] Raster file from ascii file and flattening Africa ....:)

[Corrado]

Thanks a lot Kamran!

I will test your approach and let you know ....

On Wednesday 03 December 2008 10:26:14 Kamran Safi wrote:
> Hi Corrado,
>
> Being a advanced newbie myself, I first of all understand what you mean
> by that and secondly ask you to qualify my answer.
>
> I would, tackling your problem, create a raster polygon in a metric
> equal area projection, such as Mollweide. Then you use overlay() and get
> for each polygon the set of points that are within each raster polygon.
> You need to import the xyz file in R and convert it into a Spatialpoints
> data frame.
>
> Here's the first bit
> This reads a coast line shapefile and extracts africa from it. Then uses
> the boundings boxes to produce a grid at the extent of africa. Then it
> projects that raster back to longlat for the overlay() procedure.
>
> map <- readShapePoly("E:/Science/continent.shp", ID="CONTINENT",
> proj4string = CRS("+proj=longlat"))
> africa <- as.SpatialPolygons.PolygonsList(map at polygons[1])
> africa.proj <- spTransform(africa, CRS("+proj=moll"))
> grd <- GridTopology(c(bbox(africa.proj)[1,1]+5000,
> bbox(africa.proj)[2,1]+50000), c(100000,100000),
> c(ceiling((bbox(africa.proj)[1,2] - bbox(africa.proj)[1,1]) /
> 100000),ceiling((bbox(africa.proj)[1,2] - bbox(africa.proj)[1,1]) /
> 100000)))
>
> # if you should not have a coastline of africa:
> # these are the values you'll need to produce the raster you need to
> proceed
> # bbox(africa.proj)
> #         min     max
> # r1 -2472164 6131319
> # r2 -4202811 4490010
>
>
> polys.proj <- as.SpatialPolygons.GridTopology(grd)
> proj4string(polys.proj) <- CRS("+proj=moll")
> polys <- spTransform(polys.proj, CRS("+proj=longlat"))
> # now you have a spatialpolygon in longlat proj that has equal area and
> a
> # cell size of 100km2
> #prepare a list for you results
> results <- rep(NA, length(polys at polygons))
> # use something like this to calculate the values per raster grid
> # this here is not working probably, as I didn't think about it all too
> much
> # I have though somewhere a code lying around doing exactly this step,
> so if
> # you don't succeed, let me know and I dig that out
> for(i in 1:length(polys at polygons))
> {
> results[i] <- mean(<your spatial
> points>$values[which(overlay(as.SpatialPolygons.PolygonsList(map.Proj at po
> lygons[i]), <your Spatialpoints>) != NA, ])
> }
>
> Apart from the final bits, I tested the start and it worked. The last
> bit should not be too difficult to solve. The whole thing could be more
> elegant by excluding those polygons that are in the sea. But I think
> that is something others can try to get their heads around it. Shouldn't
> be too difficult: get the centroid coordinates, overlay it with the
> costlines of Africa and convert it back into a grid...
>
> Hope this helped.
>
> Kamran
>
>
> ------------------------
> Kamran Safi
>
> Postdoctoral Research Fellow
> Institute of Zoology
> Zoological Society of London
> Regent's Park
> London NW1 4RY
>
> http://www.zoo.cam.ac.uk/ioz/people/safi.htm
>
> http://spatialr.googlepages.com
> http://asapi.wetpaint.com
>
> -----Original Message-----
> From: r-sig-geo-bounces at stat.math.ethz.ch
> [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Corrado
> Sent: 03 December 2008 09:29
> To: r-sig-ecology at r-project.org; r-sig-geo at stat.math.ethz.ch
> Subject: [R-sig-Geo] Raster file from ascii file and flattening Africa
> ....:)
>
> Dear friends,
>
> I am a kind of advanced newbie, if that makes sense.
>
> I have a text file of the form
>
> coordinate x,coordinate y,cat={real number between 250 and 450}
>
> where coordinate are expressed in latitude and longitude. The files
> represents
> measurements of the size of a skulls on sites all over Africa.
>
> >From it, I would like to build a raster file, 100 km by 100km.  There
>
> are 2
> problems:
>
> 1) Unfortunately,  in some 100km x 100km squares, there is one of the
> points
> whilst in others there are maybe 20. How do I average, so that in each
> square
> I only have 1 value representing the average?
>
> 2) How do we "flatten" Africa so that we may use 100km x 100km squares
> instead
> of 1 degree x 1 degree, without committing a geographical crime? What we
> need
> is to respect the areas ....
>
> Best regards and apologies for the silliness of the questions.



-- 
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529 at york.ac.uk