[GRASS-user] trouble with r.in.gdal to georeference a NetCDF file

Lee Eddington lee.w.eddington at gmail.com
Mon Nov 25 21:07:46 PST 2013 

I got a python script from someone on the gdal-dev list that reads in the
WRF NetCDF data and outputs a GeoTIFF file in the native Mercator
projection and I'm able to read it into GRASS OK.  But the coordinate units
are meters and I'd prefer to have degrees lat/lon.  I reposted asking about
this on gdal-dev, but have not heard back yet, and am wondering if anyone
here knows a simple way to change to angular units of degrees instead of
meters.  I've been hunting through the GRASS manual, but haven't found
anything that works except reprojecting to a lat/lon location, but then I
loose the native mercator aspect ratio of the location.  There is also a
SetAngularUnits gdal Python call that I found, but I can't get it to work.
 Here is the Python code used to transform from WRF NetCDF to GetTIFF:

from osgeo import gdal
from osgeo import osr
import numpy
import numpy.ma as ma
datafile = 'wrfout_d03_2013-11-20_12:00:00.nc'
proj_out = osr.SpatialReference()
proj_out.SetMercator(0.0, 115.02, 0.98931892612652, 0.0, 0.0)
ds_in = gdal.Open(datafile)
subdatasets = ds_in.GetSubDatasets()
variables = []
for subdataset in subdatasets:
     variables.append(subdataset[1].split(" ")[1])
     variables.append(subdataset[1].split(" ")[1])
ds_lon = gdal.Open('NETCDF:"wrfout_d03_2013-11-20_12:00:00.nc":XLONG')
ds_lat = gdal.Open('NETCDF:"wrfout_d03_2013-11-20_12:00:00.nc":XLAT')
longs = ds_lon.GetRasterBand(1).ReadAsArray()
lats = ds_lat.GetRasterBand(1).ReadAsArray()
ds_lon = None
ds_lat = None
proj_gcp = osr.SpatialReference()
proj_gcp.ImportFromEPSG(4326)
transf = osr.CoordinateTransformation(proj_gcp, proj_out)
ul = transf.TransformPoint(float(longs[0][0]), float(lats[0][0]))
lr = transf.TransformPoint(float(longs[len(longs)-1][len(longs[0])-1]),
float(lats[len(longs)-1][len(longs[0])-1]))
ur = transf.TransformPoint(float(longs[0][len(longs[0])-1]),
float(lats[0][len(longs[0])-1]))
ll = transf.TransformPoint(float(longs[len(longs)-1][0]),
float(lats[len(longs)-1][0]))
gt0 = ul[0]
gt1 = (ur[0] - gt0) / len(longs[0])
gt2 = (lr[0] - gt0 - gt1*len(longs[0])) / len(longs)
gt3 = ul[1]
gt5 = (ll[1] - gt3) / len(longs)
gt4 = (lr[1] - gt3 - gt5*len(longs) ) / len(longs[0])
gt = (gt0,gt1,gt2,gt3,gt4,gt5)
ds_u10 = gdal.Open('NETCDF:"'+datafile+'":U10')
num_bands = ds_u10.RasterCount
x_size = ds_u10.RasterXSize
y_size = ds_u10.RasterYSize
data_u10 = numpy.zeros(((num_bands, y_size, x_size)))
data_u10[:,:] = ds_u10.GetRasterBand(1).ReadAsArray()
driver = gdal.GetDriverByName( 'GTiff' )
ds_out = driver.Create( 'u10.tiff', x_size, y_size, 1, gdal.GDT_Float32 )
ds_out.SetGeoTransform( gt )
ds_out.SetProjection(proj_out.ExportToWkt())
ds_out.GetRasterBand(1).WriteArray( data_u10 )



On Sat, Nov 23, 2013 at 9:14 PM, Lee Eddington <lee.w.eddington at gmail.com>wrote:

> Markus,
>
> I can read a raster from the file into GRASS using:
>
> r.in.gdal input="NETCDF:test.nc:V10" output="V10_test3"
>
> and I get the following:
>
> GRASS 6.4.3 (bali_d03):~ > r.info V10_test3
>
>  +----------------------------------------------------------------------------+
>  | Layer:    V10_test3                      Date: Sat Nov 23 20:56:28 2013
>    |
>  | Mapset:   PERMANENT                      Login of Creator: Lee
>     |
>  | Location: bali_d03
>     |
>  | DataBase: /Users/Lee/grassdata
>     |
>  | Title:     ( V10_test3 )
>     |
>  | Timestamp: none
>    |
>
>  |----------------------------------------------------------------------------|
>  |
>    |
>  |   Type of Map:  raster               Number of Categories: 255
>     |
>  |   Data Type:    FCELL
>    |
>  |   Rows:         153
>    |
>  |   Columns:      198
>    |
>  |   Total Cells:  30294
>    |
>  |        Projection: x,y
>     |
>  |            N:        153    S:          0   Res:     1
>     |
>  |            E:        198    W:          0   Res:     1
>     |
>  |   Range of data:    min = -12.64637  max = 7.255654
>    |
>  |
>    |
>  |   Data Description:
>    |
>  |    generated by r.in.gdal
>    |
>  |
>    |
>  |   Comments:
>    |
>  |    r.in.gdal input="NETCDF:test.nc:V10" output="V10_test3"
>     |
>  |
>    |
>
>  +----------------------------------------------------------------------------+
>
> gdalinfo on test.nc gives the following (I’ve stripped out a lot of the
> metadata, but left the info the describes the Mercator projection):
>
> lees-mbp:bali Lee$ gdalinfo test.nc
> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
> Driver: netCDF/Network Common Data Format
> Files: test.nc
> Size is 512, 512
> Coordinate System is `'
> Metadata:
> .
> .
> .
>   NC_GLOBAL#CEN_LAT=-8.4095154
>   NC_GLOBAL#CEN_LON=115.02
> .
> .
> .
>   NC_GLOBAL#MAP_PROJ=3
> .
> .
> .
>   NC_GLOBAL#STAND_LON=115.02
> .
> .
>   NC_GLOBAL#TRUELAT1=-8.4095001
>   NC_GLOBAL#TRUELAT2=0
> .
> .
> .
> Subdatasets:
>   SUBDATASET_1_NAME=NETCDF:"test.nc":U10
>   SUBDATASET_1_DESC=[1x153x198] U10 (32-bit floating-point)
>   SUBDATASET_2_NAME=NETCDF:"test.nc":V10
>   SUBDATASET_2_DESC=[1x153x198] V10 (32-bit floating-point)
> Corner Coordinates:
> Upper Left  (    0.0,    0.0)
> Lower Left  (    0.0,  512.0)
> Upper Right (  512.0,    0.0)
> Lower Right (  512.0,  512.0)
> Center      (  256.0,  256.0)
>
> I can display the data in GRASS (see attachment).
>
> I have the latitude and longitude values for every raster point from
> another file and can read them into collocated rasters.  I was able to
> interactively georectify the raster, but my goal is produce automated
> scripts and it doesn’t appear that there’s a way to automate the georectify
> process (please correct me if I’m wrong).  I tried using g.setproj to
> change the projection as I know what it is, but it wouldn’t allow me to
> because the current projection is x,y.
>
> Thanks for sticking with me on this,
>
> Lee
>
>
> On Nov 23, 2013, at 7:43 AM, Markus Neteler <neteler at osgeo.org> wrote:
>
> Lee,
>
> I am not much familiar with non-map netCDF files. I do not know how to
> deal with this specific file...
> Let's hope that an answer from the community arrives. To me it doesn't
> look like a file containing maps.
>
> Best
> Markus
>
>
> On Sat, Nov 23, 2013 at 9:58 AM, Lee Eddington
> <lee.w.eddington at gmail.com> wrote:
>
> Markus,
>
> I tried and got the following:
>
> lees-mbp:bali Lee$ gdalwarp test.nc test.tif
>
> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
>
> Input file test.nc has no raster bands.
>
>
> I've attached the file in case you'd be willing to see what you can do with
> it.
>
> Thanks,
> Lee
>
>
>
> On Sun, Nov 17, 2013 at 1:29 AM, Markus Neteler <neteler at osgeo.org> wrote:
>
>
> On Wed, Nov 13, 2013 at 11:58 PM, Lee Eddington
> <lee.w.eddington at gmail.com> wrote:
>
> Markus,
>
> Could you provide me an example of using gdalwarp converting a NetCDF
> field
> to a GeoTIFF?
>
>
> I often due multiple runs here. Initially
>
> gdalwarp thefile.nc thefile.tif
>
> then I check if the target resolution is as expected (if not, see the
> -tr parameter) and so on.
> Here is an older blog page which may be of interest:
>
> http://blog.neteler.org/gdal-raster-data-tips-and-tricks/
>
> Markus
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20131125/87b11525/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: V10_test.png
Type: image/png
Size: 93806 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20131125/87b11525/attachment-0001.png>

More information about the grass-user mailing list