[GRASS-user] import of netCDF file into GRASS 7.02
Thomas Adams
tea3rd at gmail.com
Tue Jan 26 10:29:13 PST 2016
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 at 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 at 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)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20160126/12f1bc48/attachment.html>
More information about the grass-user
mailing list