<div dir="ltr">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:<div>
<br></div><div><div>from osgeo import gdal</div><div>from osgeo import osr</div><div>import numpy</div><div>import <a href="http://numpy.ma">numpy.ma</a> as ma</div><div>datafile = 'wrfout_d03_2013-11-20_12:00:<a href="http://00.nc">00.nc</a>'</div>
<div>proj_out = osr.SpatialReference()</div><div>proj_out.SetMercator(0.0, 115.02, 0.98931892612652, 0.0, 0.0)</div><div>ds_in = gdal.Open(datafile)</div><div>subdatasets = ds_in.GetSubDatasets()</div><div>variables = []</div>
<div>for subdataset in subdatasets:</div><div> variables.append(subdataset[1].split(" ")[1])</div><div> variables.append(subdataset[1].split(" ")[1])</div><div>ds_lon = gdal.Open('NETCDF:"wrfout_d03_2013-11-20_12:00:<a href="http://00.nc">00.nc</a>":XLONG')</div>
<div>ds_lat = gdal.Open('NETCDF:"wrfout_d03_2013-11-20_12:00:<a href="http://00.nc">00.nc</a>":XLAT')</div><div>longs = ds_lon.GetRasterBand(1).ReadAsArray()</div><div>lats = ds_lat.GetRasterBand(1).ReadAsArray()</div>
<div>ds_lon = None</div><div>ds_lat = None</div><div>proj_gcp = osr.SpatialReference()</div><div>proj_gcp.ImportFromEPSG(4326)</div><div>transf = osr.CoordinateTransformation(proj_gcp, proj_out)</div><div>ul = transf.TransformPoint(float(longs[0][0]), float(lats[0][0]))</div>
<div>lr = transf.TransformPoint(float(longs[len(longs)-1][len(longs[0])-1]), float(lats[len(longs)-1][len(longs[0])-1]))</div><div>ur = transf.TransformPoint(float(longs[0][len(longs[0])-1]), float(lats[0][len(longs[0])-1]))</div>
<div>ll = transf.TransformPoint(float(longs[len(longs)-1][0]), float(lats[len(longs)-1][0]))</div><div>gt0 = ul[0]</div><div>gt1 = (ur[0] - gt0) / len(longs[0])</div><div>gt2 = (lr[0] - gt0 - gt1*len(longs[0])) / len(longs)</div>
<div>gt3 = ul[1]</div><div>gt5 = (ll[1] - gt3) / len(longs)</div><div>gt4 = (lr[1] - gt3 - gt5*len(longs) ) / len(longs[0])</div><div>gt = (gt0,gt1,gt2,gt3,gt4,gt5)</div><div>ds_u10 = gdal.Open('NETCDF:"'+datafile+'":U10')</div>
<div>num_bands = ds_u10.RasterCount</div><div>x_size = ds_u10.RasterXSize</div><div>y_size = ds_u10.RasterYSize</div><div>data_u10 = numpy.zeros(((num_bands, y_size, x_size)))</div><div>data_u10[:,:] = ds_u10.GetRasterBand(1).ReadAsArray()</div>
<div>driver = gdal.GetDriverByName( 'GTiff' )</div><div>ds_out = driver.Create( 'u10.tiff', x_size, y_size, 1, gdal.GDT_Float32 )</div><div>ds_out.SetGeoTransform( gt )</div><div>ds_out.SetProjection(proj_out.ExportToWkt())</div>
<div>ds_out.GetRasterBand(1).WriteArray( data_u10 )</div></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sat, Nov 23, 2013 at 9:14 PM, Lee Eddington <span dir="ltr"><<a href="mailto:lee.w.eddington@gmail.com" target="_blank">lee.w.eddington@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">Markus,<br><br>I can read a raster from the file into GRASS using:<br><br>r.in.gdal input="NETCDF:test.nc:V10" output="V10_test3"<br>
<br>and I get the following:<br><br><font face="Courier New"><div>GRASS 6.4.3 (bali_d03):~ > <a href="http://r.info" target="_blank">r.info</a> V10_test3</div><div> +----------------------------------------------------------------------------+</div>
<div> | Layer: V10_test3 Date: Sat Nov 23 20:56:28 2013 |</div><div> | Mapset: PERMANENT Login of Creator: Lee |</div><div> | Location: bali_d03 |</div>
<div> | DataBase: /Users/Lee/grassdata |</div><div> | Title: ( V10_test3 ) |</div><div> | Timestamp: none |</div>
<div> |----------------------------------------------------------------------------|</div><div> | |</div><div> | Type of Map: raster Number of Categories: 255 |</div>
<div> | Data Type: FCELL |</div><div> | Rows: 153 |</div><div> | Columns: 198 |</div>
<div> | Total Cells: 30294 |</div><div> | Projection: x,y |</div><div> | N: 153 S: 0 Res: 1 |</div>
<div> | E: 198 W: 0 Res: 1 |</div><div> | Range of data: min = -12.64637 max = 7.255654 |</div><div> | |</div>
<div> | Data Description: |</div><div> | generated by r.in.gdal |</div><div> | |</div>
<div> | Comments: |</div><div> | r.in.gdal input="NETCDF:test.nc:V10" output="V10_test3" |</div><div> | |</div>
<div> +----------------------------------------------------------------------------+</div><div><br></div><div>gdalinfo on <a href="http://test.nc" target="_blank">test.nc</a> gives the following (I’ve stripped out a lot of the metadata, but left the info the describes the Mercator projection):</div>
<div><br></div><div><div>lees-mbp:bali Lee$ gdalinfo <a href="http://test.nc" target="_blank">test.nc</a> </div><div class="im"><div>Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute</div></div><div class="im"><div>Driver: netCDF/Network Common Data Format</div>
</div><div>Files: <a href="http://test.nc" target="_blank">test.nc</a></div><div class="im"><div>Size is 512, 512</div><div>Coordinate System is `'</div><div>Metadata:</div></div><div>.</div><div>.</div><div>.</div><div>
NC_GLOBAL#CEN_LAT=-8.4095154</div><div> NC_GLOBAL#CEN_LON=115.02</div><div>.</div></div><div>.</div><div><div>.</div><div> NC_GLOBAL#MAP_PROJ=3</div><div>.</div><div>.</div><div>.</div><div> NC_GLOBAL#STAND_LON=115.02</div>
<div>.</div><div>.</div><div> NC_GLOBAL#TRUELAT1=-8.4095001</div><div> NC_GLOBAL#TRUELAT2=0</div><div>.</div><div>.</div><div>.</div><div>Subdatasets:</div><div> SUBDATASET_1_NAME=NETCDF:"<a href="http://test.nc" target="_blank">test.nc</a>":U10</div>
<div> SUBDATASET_1_DESC=[1x153x198] U10 (32-bit floating-point)</div><div> SUBDATASET_2_NAME=NETCDF:"<a href="http://test.nc" target="_blank">test.nc</a>":V10</div><div> SUBDATASET_2_DESC=[1x153x198] V10 (32-bit floating-point)</div>
<div class="im"><div>Corner Coordinates:</div><div>Upper Left ( 0.0, 0.0)</div><div>Lower Left ( 0.0, 512.0)</div><div>Upper Right ( 512.0, 0.0)</div><div>Lower Right ( 512.0, 512.0)</div><div>Center ( 256.0, 256.0)</div>
</div></div><div><br></div></font>I can display the data in GRASS (see attachment).<img height="219" width="320" src="cid:008EC921-189D-4E41-A7D7-7AA1280E0EA7@home"><div><br></div><div>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.</div>
<div><br></div><div>Thanks for sticking with me on this,</div><div><br></div><div>Lee</div><div><div class="h5"><div><br><br>On Nov 23, 2013, at 7:43 AM, Markus Neteler <<a href="mailto:neteler@osgeo.org" target="_blank">neteler@osgeo.org</a>> wrote:<br>
<br><blockquote type="cite">Lee,<br><br>I am not much familiar with non-map netCDF files. I do not know how to<br>deal with this specific file...<br>Let's hope that an answer from the community arrives. To me it doesn't<br>
look like a file containing maps.<br><br>Best<br>Markus<br><br><br>On Sat, Nov 23, 2013 at 9:58 AM, Lee Eddington<br><<a href="mailto:lee.w.eddington@gmail.com" target="_blank">lee.w.eddington@gmail.com</a>> wrote:<br>
<blockquote type="cite">Markus,<br><br>I tried and got the following:<br><br>lees-mbp:bali Lee$ gdalwarp <a href="http://test.nc" target="_blank">test.nc</a> test.tif<br><br>Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute<br>
<br>Input file <a href="http://test.nc" target="_blank">test.nc</a> has no raster bands.<br><br><br>I've attached the file in case you'd be willing to see what you can do with<br>it.<br><br>Thanks,<br>Lee<br><br><br>
<br>On Sun, Nov 17, 2013 at 1:29 AM, Markus Neteler <<a href="mailto:neteler@osgeo.org" target="_blank">neteler@osgeo.org</a>> wrote:<br><blockquote type="cite"><br>On Wed, Nov 13, 2013 at 11:58 PM, Lee Eddington<br>
<<a href="mailto:lee.w.eddington@gmail.com" target="_blank">lee.w.eddington@gmail.com</a>> wrote:<br><blockquote type="cite">Markus,<br><br>Could you provide me an example of using gdalwarp converting a NetCDF<br>field<br>
to a GeoTIFF?<br></blockquote><br>I often due multiple runs here. Initially<br><br>gdalwarp <a href="http://thefile.nc" target="_blank">thefile.nc</a> thefile.tif<br><br>then I check if the target resolution is as expected (if not, see the<br>
-tr parameter) and so on.<br>Here is an older blog page which may be of interest:<br><br><a href="http://blog.neteler.org/gdal-raster-data-tips-and-tricks/" target="_blank">http://blog.neteler.org/gdal-raster-data-tips-and-tricks/</a><br>
<br>Markus<br></blockquote><br><br></blockquote></blockquote><br></div></div></div></div></blockquote></div><br></div>