[R-sig-Geo] projectRaster() to Goode Homolosine

[Forrest Stevens]

Hi Birgit,

I apologize for not seeing the negatives in your extent parameters,
those are indeed correct.  I think the problem that you're having is
probably occurring in however you're reading in the values and
assigning them to the raster you create.  The projectRaster() function
seems to work just fine for me going from several different
projections to Interrupted Goode Homolosine.  I don't know if this
will help you or not but I seem to be able to project rasters and
vectors easily to IGH (R 3.0.1, rgdal 0.8-9, raster 2.1-30/r2743).
Please see my test code below (which creates a header for the
downloaded goge2_0g.img file so we can read its values directly with
raster from disk) below and see if it might offer any clues:


require(raster)

## Set this to wherever you downloaded and unzipped the goge2_0g.img file:
setwd("D:/tmp/test/")

ny <- 17347
nx <- 40031

## We don't need to read the values directly as we are going to create
a header to associate with this file below:
#goge_val <- readBin("goge2_0g.img", n=ny*nx, what="integer", size=1,
signed=FALSE)

goge2_0g <- raster(xmn=-20015500,xmx=20015500,ymn=-8673500,ymx=8673500,nrows=ny,ncols=nx,crs="+proj=goode
+a=6370997")
dataType(goge2_0g) <- "INT1U"

## I wish we didn't have to do this step but creating a placeholder
BIL raster seems to be the only way to write a .hdr file:
goge2_0g <- writeStart(goge2_0g, "goge2_0g.bil", format="BIL",
datatype="INT1U", overwrite=TRUE)
writeStop(goge2_0g)
hdr(goge2_0g, format="BIL")

## Unlink (remove) our temporary .bil file leaving behind the .hdr
file for loading with raster:
unlink("goge2_0g.bil")

## Now load our raster file, which should now be in the specified
projection and map info from above:
goge2_0g <- raster("goge2_0g.img")

## Specify our projection since the BIL HDR format does not include
projection information, only coords:
projection(goge2_0g) <- "+proj=goode +a=6370997"


## Now let's try to reproject some sample data to Interrupted Goode
Homolsine and overlay it:
rs <- raster(system.file("external/test.grd", package="raster"))
rs_prj <- projectRaster(rs, crs="+proj=goode +a=6370997")

plot(goge2_0g, ext=extent(rs_prj), legend=FALSE)
plot(rs_prj, add=TRUE)



Sincerely,
Forrest

On Tue, Jun 4, 2013 at 12:44 PM, Birgit Mannig
<birgit.mannig at uni-wuerzburg.de> wrote:
> Thanks for the fast answer!
>
> This is the extent as given by the USGS -- it's a global map. Actually, what
> is given is:
>
> Projection Type: Interrupted Goode Homolosine
>     Units of measure: meters
>     Pixel Size: 1000 meters
>     Radius of sphere: 6370997 m.
>     XY corner coordinates (center of pixel) in projection units (meters):
>         Lower left: (-20015000, -8673000)
>         Upper left: (-20015000, 8673000)
>         Upper right: (20015000, 8673000)
>         Lower right: (20015000, -8673000)
>
>
> The "interrupted areas" (where there is no real map/world) are filled with
> dummy values in the binary file.
> Note the difference in sign for xmn/xmx and ymn/ymx. The resulting plot
> looks fine.
>
> All the best,
> Birgit
>
>
> Zitat von Forrest Stevens <forrest at ufl.edu>:
>
>> Hi Birgit, there's something wrong with your call to raster() as the
>> xmin/xmax and ymin/ymax values should not be identical.  Could you
>> confirm what the proper extent should be for the raster you're
>> creating?
>>
>> Forrest
>>
>> On Tue, Jun 4, 2013 at 12:27 PM, Birgit Mannig
>> <birgit.mannig at uni-wuerzburg.de> wrote:
>>>
>>> Dear all,
>>>
>>> I tried to project a Raster file (orignal projection: UTM) to a file in
>>> the
>>> Goode projection.
>>> The re-projection with proj4 was obviosly not correct, neither was a
>>> re-projection from a geographical lon/lat file.
>>> Does anyone know if this is a general problem with R/proj4?
>>>
>>> Otherwise, I might have produced a faulty raster file. I used the USGS
>>> Global Ecosystem Data (http://edc2.usgs.gov/glcc/globdoc2_0.php),
>>> provided
>>> as raw binary, and had to create the Raster myself:
>>>
>>> ny <- 17347
>>> nx <- 40031
>>> goge <-
>>>
>>> raster(xmn=-20015500,xmx=20015500,ymn=-8673500,ymx=8673500,nrows=ny,ncols=nx,crs="+proj=goode
>>> +a=6370997")
>>> values(goge) <- binfromFTP
>>>
>>> Although the created Raster File looked good to me naturally the error
>>> could
>>> still occur in this step.
>>>
>>> Afterwards I just used:
>>> projectRaster(infile,goge, method="ngb")
>>>
>>>
>>> Any hints are very much appreciated!
>>>
>>> Kind regards,
>>> Birgit
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>> --
>> Forrest R. Stevens
>> Ph.D. Candidate, QSE3 IGERT Fellow
>> Department of Geography
>> Land Use and Environmental Change Institute
>> University of Florida
>> www.clas.ufl.edu/users/forrest
>>
>
>
>
> --
> Dipl.-Geogr. Birgit Mannig
>
> Universität Würzburg
> Institut für Geographie und Geologie
> -Physische Geographie-
> Zi. 120
> Am Hubland
> 97074 Würzburg
>
> Tel. +49-931-31-81319
> Fax. +49-931-31-85544
> birgit.mannig at uni-wuerzburg.de



-- 
Forrest R. Stevens
Ph.D. Candidate, QSE3 IGERT Fellow
Department of Geography
Land Use and Environmental Change Institute
University of Florida
www.clas.ufl.edu/users/forrest