[GRASS-user] import of netCDF file into GRASS 7.02

I’m attempting to import a netCDF file (filename.nc) into grass. It is comprised of daily rainfall raster data over 730 days, for a 4 km x 4 km grid. I think the resolution is 1 km so about 16 points per day.

I’ve been digging around with gdalinfo and it says the following.

When I use r.in.gdal input=/home/jamaas/research/ra1188uea/Enigma/Rainfall_data/BC_GB_daily.nc output=BC_GB_daily

I get this error

ERROR: Projection of dataset does not appear to match current location.

when I use the -o option I get

Proceeding with import of 0 raster bands…
ERROR: Selected band (1) does not exist

I think I want all of the raster bands, i.e. for all 730 days. I must be doing something simple incorrectly.

Any suggestions most welcome.

Thanks

J

···

Dr. Jim Maas

James,

Can you provide me with a link to the dataset? I have used this data a lot before… this: https://grasswiki.osgeo.org/wiki/NetCDF may be helpful.

There is also this from Micha Silver:

These geo_em* netcdf files do not have any “proper” CRS information in the headers. So you indeed get a simple X-Y matrix with cellsize 1 when you try to import directly using GDAL.

One of the ways to work around this (explained in NCAR’s user manual, chapter 4, about p. 59) is to export to an ARCInfo ASCII file, then manually edit the ASCII file header rows and change the xllcorner, yllcorner, and cellsize to the correct values.

So first:
gdal_translate -of AAIGrid NETCDF:“geo_em.d02.nc”:GREENFRAC greenfrac.asc
Now do:
ncdump geo_em.d02.nc | grep corner
:corner_lats = 36.19099f, 39.06561f, 39.08329f, 36.20801f, 36.19086f, 39.06547f, 39.08319f, 36.20791f, 36.18649f, 39.07011f, 39.0878f, 36.20351f, 36.18636f, 39.06997f, 39.0877f, 36.2034f ;
:corner_lons = -108.7585f, -108.8684f, -104.0023f, -104.0788f, -108.7641f, -108.8742f, -103.9965f, -104.0732f, -108.7584f, -108.8686f, -104.0021f, -104.0789f, -108.7639f, -108.8744f, -103.9963f, -104.0733f ;

The above “corner_lats” and “corner_lons” are the correct Lon/Lat for the corners and center for each of the four corner pixels in the domain 2. The lower left corner of the lower left pixel is at -108.7585,36.19099. You enter these into the AAI header rows, and set the cellsize to 1000, then import that altered AAIGrid into GRASS in a WGS84 location .

One caveat: this method assumes that the whole geo_em.d03 is projected in a Lon/Lat WGS84 coordinate system, which is wrong, as you know. The way I know of to get an accurate raster is to parse out the XLAT and XLONG variables, together with the GREENFRAC variable you are interested in as a CSV list of lon,lat,greenfac values, import the values as vector points, and do an interpolation.

Another possibility comes to mind, but I’ve never tried it: You should be able to convert the above xllcorner and yllcorner values to your Lambert Conformal Conic projection in advance, and use the converted values in the AAI header rows. Then import the AAI file directly into the correct LCC location. You would do this by:

  1. create a point vector of a single point (v.in.ascii) from the Long/Lat values above in an WGS84 location
  2. project that point to your LCC location (v.proj)
  3. get the LCC X-Y coordinates (v.out.ascii)
  4. Edit the header rows with the LCC xllcorner and yllcorner values
  5. Import into an LCC location (with r.in.gdal -o as you did)

Best,
Tom

···

On Tue, Jan 26, 2016 at 11:16 AM, James Maas (MED) <J.Maas@uea.ac.uk> wrote:

I’m attempting to import a netCDF file (filename.nc) into grass. It is comprised of daily rainfall raster data over 730 days, for a 4 km x 4 km grid. I think the resolution is 1 km so about 16 points per day.

I’ve been digging around with gdalinfo and it says the following.

When I use r.in.gdal input=/home/jamaas/research/ra1188uea/Enigma/Rainfall_data/BC_GB_daily.nc output=BC_GB_daily

I get this error

ERROR: Projection of dataset does not appear to match current location.

when I use the -o option I get

Proceeding with import of 0 raster bands…
ERROR: Selected band (1) does not exist

I think I want all of the raster bands, i.e. for all 730 days. I must be doing something simple incorrectly.

Any suggestions most welcome.

Thanks

J

======================================

Driver: netCDF/Network Common Data Format
Files: BC_GB_daily.nc
BC_GB_daily.nc.aux.xml
Size is 4, 4
Coordinate System is:
PROJCS[“unnamed”,
GEOGCS[“unknown”,
DATUM[“unknown”,
SPHEROID[“Spheroid”,6377563.396,299.3249646]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433]],
PROJECTION[“Transverse_Mercator”],
PARAMETER[“latitude_of_origin”,49],
PARAMETER[“central_meridian”,0],
PARAMETER[“scale_factor”,1],
PARAMETER[“false_easting”,400],
PARAMETER[“false_northing”,-100],
UNIT[“kilometre”,1000,
AUTHORITY[“EPSG”,“9036”]]]
Origin = (328.500000000000000,526.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
crs#_CoordinateAxisTypes=GeoX GeoY
crs#_CoordinateTransformType=Projection
crs#EPSG_code=EPSG:27700
crs#false_easting=400
crs#false_northing=-100
crs#grid_mapping_name=transverse_mercator
crs#inverse_flattening=299.3249646
crs#latitude_of_projection_origin=49
crs#long_name=coordinate_reference_system
crs#longitude_of_projection_origin=-2
crs#scale_factor_at_projection_origin=0.9996012717
crs#semi_major_axis=6377563.396
crs#semi_minor_axis=6356256.91

Dr. Jim Maas


grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Thomas E Adams, III
2330 Jack Warner PKWY, #334
Tuscaloosa, AL 35401

1 (513) 739-9512 (cell)

On 26/01/16 18:16, James Maas (MED) wrote:

I'm attempting to import a netCDF file (filename.nc) into grass. It is
comprised of daily rainfall raster data over 730 days, for a 4 km x 4 km
grid. I think the resolution is 1 km so about 16 points per day.

I've been digging around with gdalinfo and it says the following.

When I use r.in.gdal
input=/home/jamaas/research/ra1188uea/Enigma/Rainfall_data/BC_GB_daily.nc output=BC_GB_daily

I would start with

gdalinfo /home/jamaas/research/ra1188uea/Enigma/Rainfall_data/BC_GB_daily.nc

To see what the dataset contains exactly, especially which bands are available.

As you can see in the wiki page that Thomas linked to, netCDF has a particular syntax in gdal, as netCDF files can contain subdatasets which in return contain bands. See also [1] for details.

Moritz

[1] http://www.gdal.org/frmt_netcdf.html

On Tue, Jan 26, 2016 at 7:29 PM, Thomas Adams <tea3rd@gmail.com> wrote:
...

One of the ways to work around this (explained in NCAR's user manual,
chapter 4, about p. 59) is to export to an ARCInfo ASCII file, then manually
edit the ASCII file header rows and change the xllcorner, yllcorner, and
cellsize to the correct values.

So first:
gdal_translate -of AAIGrid NETCDF:"geo_em.d02.nc":GREENFRAC greenfrac.asc

... did you try to run "gdalwarp" instead on the channel(s)?

Markus