Anthony, Looks close, but would be better without some of the loops. Perhaps
like this, for one file: ?
library(raster)
# I do not know the coordinates, but you can fix these:
r <- raster(nrow=448, ncol=304, xmn=0, xmx=10, ymn=0, ymx=10)
# projection(r) <- '????'
filename <- "d:/nt_20080326_f17_nrt_n.bin"
con <- file(filename, 'rb')
x <- readBin(con, "raw", 300)
x <- readBin(con,"int", size=1, signed=FALSE, 150000)
close(con)
r <- setValues(r, x)
plot(r)
sp <- as(r, 'SpatialGridDataFrame')
On Wed, Oct 12, 2011 at 1:06 PM, Anthony S Fischbach <afischb...@usgs.gov>wrote:
> I am struggling to automate this.
>
> ### script to retrieve first 3 days of data from Oct 2011 ###
> require(rgdal)
> yyyy<- "2011"
> mm<- "10"
> for(dd in 1:3){
> file <- paste("nt_",yyyy,mm,sprintf("%2.2d", dd),"_f17_nrt_n.bin",
> sep="")
> data.bin <- paste("
>
> ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/
> ", file, sep="")
> #Open binary connection to the file
> bincon<- file(data.bin, 'rb')
> #Here I tried to pass bincon to the .vrt
> ice<-readGDAL("nsidcBin.vrt")
> image(ice)
> close(bincon)
> }
>
> This fails, because I am not able to pass the binary reading connection
> (bincon) to the vrt file
> I tried to simply write the .vrt file to have the binary data explicitly
> under in the SourceFilename argument
> <SourceFilename relativetoVRT="0">
>
> ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/nt_20111011_f17_nrt_n.bin
> </SourceFilename>
> but this fails with the following error.
> in .local(.Object, ...) :
> Unable to open
>
> ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/nt_20111011_f17_nrt_n.bin
> .
> No error
>
> Would I be better off to read the binary file in as a matrix and coerce it
> into a spatial pixel dataframe from there?
>
> require(sp)
> # make loop here
> yyyy<- "2011"
> mm<- "10"
> for(dd in 1:3){
> file <- paste("nt_",yyyy,mm,sprintf("%2.2d", dd),"_f17_nrt_n.bin",
> sep="")
> data.bin <- paste("
>
> ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/
> ", file, sep="")
> bincon<- file(data.bin, 'rb')
> result <- matrix(0,448,304) #built result matrix with 448 rows
> and 304 columns
> # step past the first 300 bytes
> for(i in 1:300){readBin(bincon,"raw",1)}
> # step through the 448 rows of 304 columns
> for(i in 1:448){
> for(j in 1:304){
> x <- as.integer(readBin(bincon,"raw",1))
> if(x<251) result[i,j]<- x/250
> }
> }
> close(bincon)
> image(result)
> }
>
> >From here I would need to build the spatial pixel data frame from the
> matrix...
>
>
> Tony Fischbach, Wildlife Biologist
> Walrus Research Program
> Alaska Science Center
> U.S. Geological Survey
> 4210 University Drive
> Anchorage, AK 99508-4650
>
> Phone: 907 786-7145
> FAX: 907 786-7150
>
> afischb...@usgs.gov
> http://alaska.usgs.gov/science/biology/walrus
>
>
>
> From:
> Barry Rowlingson <b.rowling...@lancaster.ac.uk>
> To:
> Michael Sumner <mdsum...@gmail.com>
> Cc:
> Anthony Fischbach <afischb...@usgs.gov>, r-sig-geo@r-project.org
> Date:
> 10/12/2011 10:54 AM
> Subject:
> Re: [R-sig-Geo] Reading National Snow and Ice Data Center binary files
> Sent by:
> b.rowling...@googlemail.com
>
>
>
> On Wed, Oct 12, 2011 at 11:08 AM, Barry Rowlingson
> <b.rowling...@lancaster.ac.uk> wrote:
>
> > Annoyingly I don't think gdal's raw data vrt system can find out the
> > number of pixels across/down the raster - you have to code them into
> > the .vrt file. I reckon with a bit of unix (or in extremis, C code)
> > craftiness it would be possible to get the info from the header of the
> > .bin file and write the relevant .vrt file...
>
> The data format details tells is the files only come in a limited
> number of sizes. I grabbed a .bin file and wrote this .vrt file to sit
> next to it:
>
> <VRTDataset rasterXSize="304" rasterYSize="448">
> <VRTRasterBand dataType="Byte" band="1" subClass="VRTRawRasterBand">
> <SourceFilename
> relativetoVRT="1">nt_20111010_f17_nrt_n.bin</SourceFilename>
> <ImageOffset>300</ImageOffset>
> <PixelOffset>1</PixelOffset>
> <LineOffset>304</LineOffset>
> </VRTRasterBand>
> </VRTDataset>
>
> Then in R:
>
> require(rgdal)
> ice=readGDAL("nt_20111010_f17_nrt_n.vrt")
> image(ice)
>
> and there's Greenland, and various other bits of rapidly melting ice...
>
> Hope this helps - the important bits of the .vrt file are the 300
> offset the raster X and Y size and LineOffset(which =XSize) and
> obviously the filename...
>
> I'd write a script to create these .vrt files and then convert
> everything to a standard format like GeoTiff... Note you'll have to
> add the projection info on manually... etc etc etc.
>
> Barry
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
