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

Reply via email to