read.topo() should handle NOAA netcdf files

[dankelley]

This came up in my testing of download.topo() ... the ascii files are really large, and I don't see the real point in encouraging people to use them, since they are likely to be reading them with read.topo() anyway.

[dankelley]

Hm. Well, it's supposed to! But I see that the files served up from http://maps.ngdc.noaa.gov/viewers/wcs-client/ have different variable names to those that read.topo() is expected. I'll need to sniff around some files on my computer, as well, so I can revise read.topo() to handle these two formats. I'll make a comment for each type, following this.

[dankelley]

The new type (server http://maps.ngdc.noaa.gov/viewers/wcs-client/ with format=netcdf) has the following ncdf content. The present read.topo() code is looking for x_range, y_range and z ... clearly not here! I'm going to try the other netcdf variant on the server next.

 ncdf
File ~/data/topo/topo_60W_40E_50S_50N_4min.nc (NC_FORMAT_CLASSIC):

     2 variables (excluding dimension variables):
        char crs[]   
            grid_mapping_name: latitude_longitude
            longitude_of_prime_meridian: 0
            semi_major_axis: 6378137
            inverse_flattening: 298.257223563
            spatial_ref: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
            GeoTransform: -60.03363350000001 0.066667 0 50.03358350000001 0 -0.066667 
        short Band1[lon,lat]   
            long_name: GDAL Band Number 1
            _FillValue: -32767
            grid_mapping: crs

     2 dimensions:
        lon  Size:1501
            standard_name: longitude
            long_name: longitude
            units: degrees_east
        lat  Size:1501
            standard_name: latitude
            long_name: latitude
            units: degrees_north

    6 global attributes:
        GDAL_TIFFTAG_RESOLUTIONUNIT: 2
        GDAL_TIFFTAG_XRESOLUTION: 72
        GDAL_TIFFTAG_YRESOLUTION: 72
        Conventions: CF-1.5
        GDAL: GDAL 1.11.0, released 2014/04/16
        history: Fri Aug 19 12:05:00 2016: GDAL CreateCopy( /nfs/www_tmp/mapserver/57b74a4c_17d1_0.nc, ... )

[dankelley]

The server has another type labelled GMT NetCDF in the pulldown menu. It codes to format=gmt in the URL. And, oh, I see now, this format has the expected entries in the netcdf file. I am going to disallow the "netcdf" format from download.topo() for now, to avoid having to recode things. (I plan to look into the data files to see if the results are different, though.)

[dankelley]

Oh, lovely. The geometry is coming out wrong.

Error in as.topo(longitude, latitude, z, filename = file) : 
  longitude vector has length 1502, which does not match matrix width 1501

PS. these babblings are just notes to myself -- I hope other devels will ignore them!!

[dankelley]

I'm guessing maybe the gmt format is for a centred grid.

[dankelley]

The file I downloaded in the other format (the old format) has as below

File ~/data/topo/topo_60W_40E_50S_50N_4min_gmt.nc (NC_FORMAT_CLASSIC) (GMT format):

     6 variables (excluding dimension variables):
        double x_range[side]   
            units: meters
        double y_range[side]   
            units: meters
        double z_range[side]   
            units: meters
        double spacing[side]   
        int dimension[side]   
        short z[xysize]   
            scale_factor: 1
            add_offset: 0
            node_offset: 1

     2 dimensions:
        side  Size:2
        xysize  Size:2253001

    2 global attributes:
        title: 
        source: 

[dankelley]

I have quite a lot of the code working, and am hoping to reformulate the topoWorld dataset with a download using this. (At the moment, this dataset is created by a hacked-together script that reads a binary file.) I am noticing some odd things with the NOAA website, though.

1. It permits east=-180 and west=180 and it gives different results for these two locations, even though they are the same spot on the earth.
2. I cannot find any documentation on whether the data that are served up by NOAA are cell-centred or grid-centred. I emailed the support team today, to enquire about that. I would like to be quite clear on this, especially for decimated data (i.e. would the offset by 1/2 arc-second as in the original file or 1/2 of whatever grid spacing we are requesting??)

[dankelley]

I'm closing this, since it's been 4 days since I asked tech support my question on how the grids are defined. If they answer me, I'll reopen this and continue on. For now, though, there is nothing I feel inclined to do, lacking documentation.