<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Nov 26, 2013 at 8:32 AM, 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 dir="ltr">Roger,<div><br></div><div>Attached is the script I'm using.</div><div><br></div><div>Here are the results from gdalinfo on the output file:</div>
<div><br></div><div><div>lees-mbp:full Lee$ gdalinfo u10_2013112018.tif</div>
<div>Driver: GTiff/GeoTIFF</div><div>Files: u10_2013112018.tif</div><div>Size is 198, 153</div><div>Coordinate System is:</div><div>PROJCS["unnamed",</div><div> GEOGCS["unknown",</div><div> DATUM["unknown",</div>
<div> SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],</div><div> PRIMEM["Greenwich",0],</div><div> UNIT["degree",0.0174532925199433]],</div><div> PROJECTION["Mercator_1SP"],</div>
<div> PARAMETER["central_meridian",115.02],</div><div> PARAMETER["scale_factor",0.9893189261265199],</div><div> PARAMETER["false_easting",0],</div><div> PARAMETER["false_northing",0],</div>
<div> UNIT["metre",1,</div><div> AUTHORITY["EPSG","9001"]]]</div><div>Origin = (-295898.731194700172637,-696504.655801336281002)</div><div>Pixel Size = (2988.872338323183612,-2964.865850256682734)</div>
<div>Metadata:</div><div> AREA_OR_POINT=Area</div><div>Image Structure Metadata:</div><div> INTERLEAVE=BAND</div><div>Corner Coordinates:</div><div>Upper Left ( -295898.731, -696504.656) (112d19'59.51"E, 6d21'13.48"S)</div>
<div>Lower Left ( -295898.731,-<a href="tel:1150129.131" value="+551150129131" target="_blank">1150129.131</a>) (112d19'59.51"E, 10d27'15.97"S)</div><div>Upper Right ( 295897.992, -696504.656) (117d42'24.46"E, 6d21'13.48"S)</div>
<div>Lower Right ( 295897.992,-<a href="tel:1150129.131" value="+551150129131" target="_blank">1150129.131</a>) (117d42'24.46"E, 10d27'15.97"S)</div>
<div>Center ( -0.370, -923316.893) (115d 1'11.99"E, 8d24'34.51"S)</div><div>Band 1 Block=198x10 Type=Float32, ColorInterp=Gray</div></div><div><br></div><div><br></div><div>When I import this into GRASS I get coordinates of meters instead of lat/lon degrees. I want to keep the mercator projection of the model and not convert to a lat/lon projection.</div>
</div></blockquote><div><br></div><div>Mercator projection is by definition in linear units (meters, km, whatever), not angular units (degrees). </div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">
<div><br></div><div>Lee</div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Nov 26, 2013 at 1:45 AM, Roger Veciana i Rovira <span dir="ltr"><<a href="mailto:rveciana@gmail.com" target="_blank">rveciana@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 dir="ltr">Lee, <div><br></div><div>If I understand what are you saiying, you have to define a proper projection first:</div>
<div><br></div><div>proj_out = osr.SpatialReference()<br></div><div>proj_out.ImportFromEPSG(4326)</div>
<div><br></div><div>Which will create the lat lon projection you need (using WGS84 datum, maybe this won't be exact in your model). Then, you have to set the geotransform directly in degrees as you want.</div><span><font color="#888888"><div>
<br>
</div><div>Roger</div></font></span></div><div><div class="h5"><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">2013/11/26 Lee Eddington <span dir="ltr"><<a href="mailto:lee.w.eddington@gmail.com" target="_blank">lee.w.eddington@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">I think I might be getting close. I found the following which I'm hoping will create the GeoTiff with angular units of degrees instead of linear units of meters:<div>
<br></div>OGRErr OGRSpatialReference::SetAngularUnits(const char * pszUnitsName,<br>
double dfInRadians <br>)<br><br><div>and tried to use it, but don't know the proper pszUnitsName. I've tried "degree", "degrees", "DEGREE", "DEGREES", "DEG", "SRS_UA_DEGREE", but always get a return code of 6 instead of 0. I can't find a list of names or any examples of how to use this properly. </div>
</div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Nov 25, 2013 at 2:00 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 dir="ltr">Roger,<div><br></div><div>Thank you for the info and link. Using a modified version of your code I was able to create and import into GRASS a GeoTIFF with the proper georeferencing. The coordinate system units are in meters though, and I'd prefer to have them in degrees lat/lon. Do you know if there's a way to specify that using the oar or gdal routines in your code?</div>
<div><br></div><div>Thanks,</div><div>Lee</div><div><br></div></div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Nov 25, 2013 at 10:23 AM, Etienne Tourigny <span dir="ltr"><<a href="mailto:etourigny.dev@gmail.com" target="_blank">etourigny.dev@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 dir="ltr">Forwarding to the list and others...<br><br><div class="gmail_quote">---------- Forwarded message ----------<br>
From: <b class="gmail_sendername">Roger Veciana i Rovira</b> <span dir="ltr"><<a href="mailto:rveciana@gmail.com" target="_blank">rveciana@gmail.com</a>></span><br>
Date: Mon, Nov 25, 2013 at 9:35 AM<br>Subject: Re: [gdal-dev] WRF NetCDF import<br>To: Etienne Tourigny <<a href="mailto:etourigny.dev@gmail.com" target="_blank">etourigny.dev@gmail.com</a>><br><br><br><div dir="ltr">
As Etienne told you, the geotransform and projection from your files is not in the format that GDAL can understand.<div>
I wrote a script that reads the projection and calculates the geotransform to save it in a GeoTIFF file:</div>
<div><br></div><div><a href="http://geoexamples.blogspot.com.es/2013/09/reading-wrf-netcdf-files-with-gdal.html" target="_blank">http://geoexamples.blogspot.com.es/2013/09/reading-wrf-netcdf-files-with-gdal.html</a><br></div>
<div><br></div>
<div>The script also interpolates the WRF sigma levels to the desired pressure level, so you will have to look for the functions that you need.</div><div><br></div><div>Anyway, the solution is quite complicated, and using cdo, nco or other tools (I use ARWPost) may be easier.</div>
<span><font color="#888888">
<div><br></div><div>Roger</div></font></span></div><div><div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">2013/11/24 Etienne Tourigny <span dir="ltr"><<a href="mailto:etourigny.dev@gmail.com" target="_blank">etourigny.dev@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>The "<span style="font-family:arial,sans-serif;font-size:13px">georeferencing info</span>" you refer to is non-standard, and GDAL cannot use that. It (and presumably the grass provider, although I don't know) use the CF standard for encoding lat/lon information. what is the output of </div>
<div>gdalinfo NETCDF:test_full.nc:HGT</div><div><br></div><div>I don't know how the grass provider deals with subdatasets, but if you open the file in qgis (using the gdal provider, it will sllow you to select the subdataset you want).</div>
<div><br></div><div>You could also extract the variable with another tool like cdo or nco.</div><div><br></div><div>Etienne</div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote"><div><div>
On Sat, Nov 23, 2013 at 7:47 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>
</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div><div dir="ltr">I'm trying to use r.in.gdal in GRASS GIS to import Weather Research & Forecasting (WRF) model NetCDF output. There have been some tips provided by the GRASS users list, but none have worked for me. I can import the data as a simple x,y array with no georeferencing, but information about georeferencing is in the file as other programs read, georeference and display the file correctly.<div>
<br></div><div>gdalinfo produces the following:</div><div><br></div><div>lees-mbp:full Lee$ gdalinfo <a href="http://test_full.nc" target="_blank">test_full.nc</a></div><div>Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute</div>
<div>
Driver: netCDF/Network Common Data Format</div><div>Files: <a href="http://test_full.nc" target="_blank">test_full.nc</a></div><div>Size is 512, 512</div><div>Coordinate System is `'</div><div>Metadata:</div><div> NC_GLOBAL#BL_PBL_PHYSICS=1</div>
<div> NC_GLOBAL#BOTTOM-TOP_GRID_DIMENSION=28</div><div> NC_GLOBAL#BOTTOM-TOP_PATCH_END_STAG=28</div><div> NC_GLOBAL#BOTTOM-TOP_PATCH_END_UNSTAG=27</div><div> NC_GLOBAL#BOTTOM-TOP_PATCH_START_STAG=1</div><div> NC_GLOBAL#BOTTOM-TOP_PATCH_START_UNSTAG=1</div>
<div> NC_GLOBAL#BUCKET_J=-1</div><div> NC_GLOBAL#BUCKET_MM=-1</div><div> NC_GLOBAL#CEN_LAT=-8.4095154</div><div> NC_GLOBAL#CEN_LON=115.02</div><div> NC_GLOBAL#CU_PHYSICS=0</div><div> NC_GLOBAL#DAMP_OPT=0</div><div>
NC_GLOBAL#DAMPCOEF=0.2</div>
<div> NC_GLOBAL#DFI_OPT=0</div><div> NC_GLOBAL#DIFF_6TH_FACTOR=0.12</div><div> NC_GLOBAL#DIFF_6TH_OPT=0</div><div> NC_GLOBAL#DIFF_OPT=1</div><div> NC_GLOBAL#DT=16.666666</div><div> NC_GLOBAL#DX=3000</div><div> NC_GLOBAL#DY=3000</div>
<div> NC_GLOBAL#FEEDBACK=0</div><div> NC_GLOBAL#GFDDA_END_H=0</div><div> NC_GLOBAL#GFDDA_INTERVAL_M=0</div><div> NC_GLOBAL#GMT=12</div><div> NC_GLOBAL#GRID_FDDA=0</div><div> NC_GLOBAL#GRID_ID=3</div><div> NC_GLOBAL#GRID_SFDDA=0</div>
<div> NC_GLOBAL#GRIDTYPE=C</div><div> NC_GLOBAL#HYPSOMETRIC_OPT=2</div><div> NC_GLOBAL#I_PARENT_START=34</div><div> NC_GLOBAL#ISFTCFLX=0</div><div> NC_GLOBAL#ISHALLOW=0</div><div> NC_GLOBAL#ISICE=24</div><div> NC_GLOBAL#ISLAKE=-1</div>
<div> NC_GLOBAL#ISOILWATER=14</div><div> NC_GLOBAL#ISURBAN=1</div><div> NC_GLOBAL#ISWATER=16</div><div> NC_GLOBAL#J_PARENT_START=34</div><div> NC_GLOBAL#JULDAY=324</div><div> NC_GLOBAL#JULYR=2013</div><div> NC_GLOBAL#KHDIF=0</div>
<div> NC_GLOBAL#KM_OPT=4</div><div> NC_GLOBAL#KVDIF=0</div><div> NC_GLOBAL#MAP_PROJ=3</div><div> NC_GLOBAL#MFSHCONV=0</div><div> NC_GLOBAL#MMINLU=USGS</div><div> NC_GLOBAL#MOAD_CEN_LAT=-8.4095078</div><div> NC_GLOBAL#MOIST_ADV_OPT=1</div>
<div> NC_GLOBAL#MP_PHYSICS=3</div><div> NC_GLOBAL#NUM_LAND_CAT=24</div><div> NC_GLOBAL#OBS_NUDGE_OPT=0</div><div> NC_GLOBAL#OMLCALL=0</div><div> NC_GLOBAL#PARENT_GRID_RATIO=3</div><div> NC_GLOBAL#PARENT_ID=2</div><div>
NC_GLOBAL#POLE_LAT=90</div><div> NC_GLOBAL#POLE_LON=0</div><div> NC_GLOBAL#PREC_ACC_DT=0</div><div> NC_GLOBAL#RA_LW_PHYSICS=1</div><div> NC_GLOBAL#RA_SW_PHYSICS=1</div><div> NC_GLOBAL#SCALAR_ADV_OPT=1</div><div> NC_GLOBAL#SF_SFCLAY_PHYSICS=1</div>
<div> NC_GLOBAL#SF_SURFACE_PHYSICS=2</div><div> NC_GLOBAL#SF_URBAN_PHYSICS=0</div><div> NC_GLOBAL#SGFDDA_END_H=0</div><div> NC_GLOBAL#SGFDDA_INTERVAL_M=0</div><div> NC_GLOBAL#SHCU_PHYSICS=0</div><div> NC_GLOBAL#SIMULATION_START_DATE=2013-11-20_12:00:00</div>
<div> NC_GLOBAL#SMOOTH_OPTION=0</div><div> NC_GLOBAL#SOUTH-NORTH_GRID_DIMENSION=154</div><div> NC_GLOBAL#SOUTH-NORTH_PATCH_END_STAG=154</div><div> NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG=153</div><div> NC_GLOBAL#SOUTH-NORTH_PATCH_START_STAG=1</div>
<div> NC_GLOBAL#SOUTH-NORTH_PATCH_START_UNSTAG=1</div><div> NC_GLOBAL#SST_UPDATE=0</div><div> NC_GLOBAL#STAND_LON=115.02</div><div> NC_GLOBAL#START_DATE=2013-11-20_12:00:00</div><div> NC_GLOBAL#SURFACE_INPUT_SOURCE=1</div>
<div> NC_GLOBAL#SWRAD_SCAT=1</div><div> NC_GLOBAL#TITLE= OUTPUT FROM WRF V3.4 MODEL</div><div> NC_GLOBAL#TKE_ADV_OPT=1</div><div> NC_GLOBAL#TRUELAT1=-8.4095001</div><div> NC_GLOBAL#TRUELAT2=0</div><div> NC_GLOBAL#W_DAMPING=0</div>
<div> NC_GLOBAL#WEST-EAST_GRID_DIMENSION=199</div><div> NC_GLOBAL#WEST-EAST_PATCH_END_STAG=199</div><div> NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG=198</div><div> NC_GLOBAL#WEST-EAST_PATCH_START_STAG=1</div><div> NC_GLOBAL#WEST-EAST_PATCH_START_UNSTAG=1</div>
<div>Subdatasets:</div><div> SUBDATASET_1_NAME=NETCDF:"<a href="http://test_full.nc" target="_blank">test_full.nc</a>":Times</div><div> SUBDATASET_1_DESC=[1x19] Times (8-bit character)</div><div> SUBDATASET_2_NAME=NETCDF:"<a href="http://test_full.nc" target="_blank">test_full.nc</a>":LU_INDEX</div>
<div> SUBDATASET_2_DESC=[1x153x198] LU_INDEX (32-bit floating-point)</div><div><br></div><div>(more SUBDATASETs)</div><div><br></div><div> SUBDATASET_106_NAME=NETCDF:"<a href="http://test_full.nc" target="_blank">test_full.nc</a>":LANDMASK</div>
<div> SUBDATASET_106_DESC=[1x153x198] LANDMASK (32-bit floating-point)</div><div> SUBDATASET_107_NAME=NETCDF:"<a href="http://test_full.nc" target="_blank">test_full.nc</a>":SST</div><div> SUBDATASET_107_DESC=[1x153x198] SST (32-bit floating-point)</div>
<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><br></div><div>My first question is why am I getting the following warning?:</div><div><br></div><div>Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute<br></div><div><br></div><div>Next, I see what appears is georeferencing info in the following lines:</div>
<div><br></div><div><div> NC_GLOBAL#CEN_LAT=-8.4095154</div><div> NC_GLOBAL#CEN_LON=115.02</div></div><div><div> NC_GLOBAL#MAP_PROJ=3 (Mercator)<br></div></div><div><div> NC_GLOBAL#MOAD_CEN_LAT=-8.4095078</div></div>
<div>
<div> NC_GLOBAL#STAND_LON=115.02</div></div><div><div> NC_GLOBAL#TRUELAT1=-8.4095001</div></div><div><br></div><div>but obviously this is not being interpreted or used.</div><div><br></div><div>Trying to use gdalwarp:</div>
<div><br></div>lees-mbp:full Lee$ gdalwarp NETCDF:test_full.nc:HGT test_full_HGT.tif<br>Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute<br>ERROR 1: Unable to compute a transformation between pixel/line<br>and georeferenced coordinates for NETCDF:test_full.nc:HGT.<br>
There is no affine transformation and no GCPs.<div><br></div><div>Can anyone explain what I need to do?</div><div><br></div><div>Thanks,</div><div>Lee</div><div><br></div></div>
<br></div></div>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br></blockquote></div><br></div>
<br>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br></blockquote></div><br></div>
</div></div></div></div></div><br></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</div></div></div></div></blockquote></div><br></div>
<br>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br></blockquote></div><br></div></div>