[gdal-dev] WRF NetCDF import

Etienne Tourigny etourigny.dev at gmail.com
Tue Nov 26 13:49:38 PST 2013 

On Tue, Nov 26, 2013 at 8:32 AM, Lee Eddington <lee.w.eddington at gmail.com>wrote:

> Roger,
>
> Attached is the script I'm using.
>
> Here are the results from gdalinfo on the output file:
>
> lees-mbp:full Lee$ gdalinfo u10_2013112018.tif
> Driver: GTiff/GeoTIFF
> Files: u10_2013112018.tif
> Size is 198, 153
> Coordinate System is:
> PROJCS["unnamed",
>     GEOGCS["unknown",
>         DATUM["unknown",
>             SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
>         PRIMEM["Greenwich",0],
>         UNIT["degree",0.0174532925199433]],
>     PROJECTION["Mercator_1SP"],
>     PARAMETER["central_meridian",115.02],
>     PARAMETER["scale_factor",0.9893189261265199],
>     PARAMETER["false_easting",0],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]]]
> Origin = (-295898.731194700172637,-696504.655801336281002)
> Pixel Size = (2988.872338323183612,-2964.865850256682734)
> Metadata:
>   AREA_OR_POINT=Area
> Image Structure Metadata:
>   INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left  ( -295898.731, -696504.656) (112d19'59.51"E,  6d21'13.48"S)
> Lower Left  ( -295898.731,-1150129.131) (112d19'59.51"E, 10d27'15.97"S)
> Upper Right (  295897.992, -696504.656) (117d42'24.46"E,  6d21'13.48"S)
> Lower Right (  295897.992,-1150129.131) (117d42'24.46"E, 10d27'15.97"S)
> Center      (      -0.370, -923316.893) (115d 1'11.99"E,  8d24'34.51"S)
> Band 1 Block=198x10 Type=Float32, ColorInterp=Gray
>
>
> 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.
>

Mercator projection is by definition in linear units (meters, km,
whatever), not angular units (degrees).


> Lee
>
>
>
> On Tue, Nov 26, 2013 at 1:45 AM, Roger Veciana i Rovira <
> rveciana at gmail.com> wrote:
>
>> Lee,
>>
>> If I understand what are you saiying, you have to define a proper
>> projection first:
>>
>> proj_out = osr.SpatialReference()
>> proj_out.ImportFromEPSG(4326)
>>
>> 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.
>>
>> Roger
>>
>>
>> 2013/11/26 Lee Eddington <lee.w.eddington at gmail.com>
>>
>>> 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:
>>>
>>> OGRErr OGRSpatialReference::SetAngularUnits(const char * pszUnitsName,
>>> double dfInRadians
>>> )
>>>
>>> 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.
>>>
>>>
>>> On Mon, Nov 25, 2013 at 2:00 PM, Lee Eddington <
>>> lee.w.eddington at gmail.com> wrote:
>>>
>>>> Roger,
>>>>
>>>> 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?
>>>>
>>>> Thanks,
>>>> Lee
>>>>
>>>>
>>>>
>>>> On Mon, Nov 25, 2013 at 10:23 AM, Etienne Tourigny <
>>>> etourigny.dev at gmail.com> wrote:
>>>>
>>>>> Forwarding to the list and others...
>>>>>
>>>>> ---------- Forwarded message ----------
>>>>> From: Roger Veciana i Rovira <rveciana at gmail.com>
>>>>> Date: Mon, Nov 25, 2013 at 9:35 AM
>>>>> Subject: Re: [gdal-dev] WRF NetCDF import
>>>>> To: Etienne Tourigny <etourigny.dev at gmail.com>
>>>>>
>>>>>
>>>>> As Etienne told you, the geotransform and projection from your files
>>>>> is not in the format that GDAL can understand.
>>>>> I wrote a script that reads the projection and calculates the
>>>>> geotransform to save it in a GeoTIFF file:
>>>>>
>>>>>
>>>>> http://geoexamples.blogspot.com.es/2013/09/reading-wrf-netcdf-files-with-gdal.html
>>>>>
>>>>> 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.
>>>>>
>>>>> Anyway, the solution is quite complicated, and using cdo, nco or other
>>>>> tools (I use ARWPost) may be easier.
>>>>>
>>>>> Roger
>>>>>
>>>>>
>>>>> 2013/11/24 Etienne Tourigny <etourigny.dev at gmail.com>
>>>>>
>>>>>> The "georeferencing info" 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
>>>>>> gdalinfo NETCDF:test_full.nc:HGT
>>>>>>
>>>>>> 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).
>>>>>>
>>>>>> You could also extract the variable with another tool like cdo or nco.
>>>>>>
>>>>>> Etienne
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Sat, Nov 23, 2013 at 7:47 PM, Lee Eddington <
>>>>>> lee.w.eddington at gmail.com> wrote:
>>>>>>
>>>>>>> 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.
>>>>>>>
>>>>>>> gdalinfo produces the following:
>>>>>>>
>>>>>>> lees-mbp:full Lee$ gdalinfo test_full.nc
>>>>>>> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
>>>>>>>  Driver: netCDF/Network Common Data Format
>>>>>>> Files: test_full.nc
>>>>>>> Size is 512, 512
>>>>>>> Coordinate System is `'
>>>>>>> Metadata:
>>>>>>>   NC_GLOBAL#BL_PBL_PHYSICS=1
>>>>>>>   NC_GLOBAL#BOTTOM-TOP_GRID_DIMENSION=28
>>>>>>>   NC_GLOBAL#BOTTOM-TOP_PATCH_END_STAG=28
>>>>>>>   NC_GLOBAL#BOTTOM-TOP_PATCH_END_UNSTAG=27
>>>>>>>   NC_GLOBAL#BOTTOM-TOP_PATCH_START_STAG=1
>>>>>>>   NC_GLOBAL#BOTTOM-TOP_PATCH_START_UNSTAG=1
>>>>>>>   NC_GLOBAL#BUCKET_J=-1
>>>>>>>   NC_GLOBAL#BUCKET_MM=-1
>>>>>>>   NC_GLOBAL#CEN_LAT=-8.4095154
>>>>>>>   NC_GLOBAL#CEN_LON=115.02
>>>>>>>   NC_GLOBAL#CU_PHYSICS=0
>>>>>>>   NC_GLOBAL#DAMP_OPT=0
>>>>>>>   NC_GLOBAL#DAMPCOEF=0.2
>>>>>>>   NC_GLOBAL#DFI_OPT=0
>>>>>>>   NC_GLOBAL#DIFF_6TH_FACTOR=0.12
>>>>>>>   NC_GLOBAL#DIFF_6TH_OPT=0
>>>>>>>   NC_GLOBAL#DIFF_OPT=1
>>>>>>>   NC_GLOBAL#DT=16.666666
>>>>>>>   NC_GLOBAL#DX=3000
>>>>>>>   NC_GLOBAL#DY=3000
>>>>>>>   NC_GLOBAL#FEEDBACK=0
>>>>>>>   NC_GLOBAL#GFDDA_END_H=0
>>>>>>>   NC_GLOBAL#GFDDA_INTERVAL_M=0
>>>>>>>   NC_GLOBAL#GMT=12
>>>>>>>   NC_GLOBAL#GRID_FDDA=0
>>>>>>>   NC_GLOBAL#GRID_ID=3
>>>>>>>   NC_GLOBAL#GRID_SFDDA=0
>>>>>>>   NC_GLOBAL#GRIDTYPE=C
>>>>>>>   NC_GLOBAL#HYPSOMETRIC_OPT=2
>>>>>>>   NC_GLOBAL#I_PARENT_START=34
>>>>>>>   NC_GLOBAL#ISFTCFLX=0
>>>>>>>   NC_GLOBAL#ISHALLOW=0
>>>>>>>   NC_GLOBAL#ISICE=24
>>>>>>>   NC_GLOBAL#ISLAKE=-1
>>>>>>>   NC_GLOBAL#ISOILWATER=14
>>>>>>>   NC_GLOBAL#ISURBAN=1
>>>>>>>   NC_GLOBAL#ISWATER=16
>>>>>>>   NC_GLOBAL#J_PARENT_START=34
>>>>>>>   NC_GLOBAL#JULDAY=324
>>>>>>>   NC_GLOBAL#JULYR=2013
>>>>>>>   NC_GLOBAL#KHDIF=0
>>>>>>>   NC_GLOBAL#KM_OPT=4
>>>>>>>   NC_GLOBAL#KVDIF=0
>>>>>>>   NC_GLOBAL#MAP_PROJ=3
>>>>>>>   NC_GLOBAL#MFSHCONV=0
>>>>>>>   NC_GLOBAL#MMINLU=USGS
>>>>>>>   NC_GLOBAL#MOAD_CEN_LAT=-8.4095078
>>>>>>>   NC_GLOBAL#MOIST_ADV_OPT=1
>>>>>>>   NC_GLOBAL#MP_PHYSICS=3
>>>>>>>   NC_GLOBAL#NUM_LAND_CAT=24
>>>>>>>   NC_GLOBAL#OBS_NUDGE_OPT=0
>>>>>>>   NC_GLOBAL#OMLCALL=0
>>>>>>>   NC_GLOBAL#PARENT_GRID_RATIO=3
>>>>>>>   NC_GLOBAL#PARENT_ID=2
>>>>>>>   NC_GLOBAL#POLE_LAT=90
>>>>>>>   NC_GLOBAL#POLE_LON=0
>>>>>>>   NC_GLOBAL#PREC_ACC_DT=0
>>>>>>>   NC_GLOBAL#RA_LW_PHYSICS=1
>>>>>>>   NC_GLOBAL#RA_SW_PHYSICS=1
>>>>>>>   NC_GLOBAL#SCALAR_ADV_OPT=1
>>>>>>>   NC_GLOBAL#SF_SFCLAY_PHYSICS=1
>>>>>>>   NC_GLOBAL#SF_SURFACE_PHYSICS=2
>>>>>>>   NC_GLOBAL#SF_URBAN_PHYSICS=0
>>>>>>>   NC_GLOBAL#SGFDDA_END_H=0
>>>>>>>   NC_GLOBAL#SGFDDA_INTERVAL_M=0
>>>>>>>   NC_GLOBAL#SHCU_PHYSICS=0
>>>>>>>   NC_GLOBAL#SIMULATION_START_DATE=2013-11-20_12:00:00
>>>>>>>   NC_GLOBAL#SMOOTH_OPTION=0
>>>>>>>   NC_GLOBAL#SOUTH-NORTH_GRID_DIMENSION=154
>>>>>>>   NC_GLOBAL#SOUTH-NORTH_PATCH_END_STAG=154
>>>>>>>   NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG=153
>>>>>>>   NC_GLOBAL#SOUTH-NORTH_PATCH_START_STAG=1
>>>>>>>   NC_GLOBAL#SOUTH-NORTH_PATCH_START_UNSTAG=1
>>>>>>>   NC_GLOBAL#SST_UPDATE=0
>>>>>>>   NC_GLOBAL#STAND_LON=115.02
>>>>>>>   NC_GLOBAL#START_DATE=2013-11-20_12:00:00
>>>>>>>   NC_GLOBAL#SURFACE_INPUT_SOURCE=1
>>>>>>>   NC_GLOBAL#SWRAD_SCAT=1
>>>>>>>   NC_GLOBAL#TITLE= OUTPUT FROM WRF V3.4 MODEL
>>>>>>>   NC_GLOBAL#TKE_ADV_OPT=1
>>>>>>>   NC_GLOBAL#TRUELAT1=-8.4095001
>>>>>>>   NC_GLOBAL#TRUELAT2=0
>>>>>>>   NC_GLOBAL#W_DAMPING=0
>>>>>>>   NC_GLOBAL#WEST-EAST_GRID_DIMENSION=199
>>>>>>>   NC_GLOBAL#WEST-EAST_PATCH_END_STAG=199
>>>>>>>   NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG=198
>>>>>>>   NC_GLOBAL#WEST-EAST_PATCH_START_STAG=1
>>>>>>>   NC_GLOBAL#WEST-EAST_PATCH_START_UNSTAG=1
>>>>>>> Subdatasets:
>>>>>>>   SUBDATASET_1_NAME=NETCDF:"test_full.nc":Times
>>>>>>>   SUBDATASET_1_DESC=[1x19] Times (8-bit character)
>>>>>>>   SUBDATASET_2_NAME=NETCDF:"test_full.nc":LU_INDEX
>>>>>>>   SUBDATASET_2_DESC=[1x153x198] LU_INDEX (32-bit floating-point)
>>>>>>>
>>>>>>> (more SUBDATASETs)
>>>>>>>
>>>>>>>   SUBDATASET_106_NAME=NETCDF:"test_full.nc":LANDMASK
>>>>>>>   SUBDATASET_106_DESC=[1x153x198] LANDMASK (32-bit floating-point)
>>>>>>>   SUBDATASET_107_NAME=NETCDF:"test_full.nc":SST
>>>>>>>   SUBDATASET_107_DESC=[1x153x198] SST (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)
>>>>>>>
>>>>>>> My first question is why am I getting the following warning?:
>>>>>>>
>>>>>>> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
>>>>>>>
>>>>>>> Next, I see what appears is georeferencing info in the following
>>>>>>> lines:
>>>>>>>
>>>>>>>   NC_GLOBAL#CEN_LAT=-8.4095154
>>>>>>>   NC_GLOBAL#CEN_LON=115.02
>>>>>>>   NC_GLOBAL#MAP_PROJ=3 (Mercator)
>>>>>>>   NC_GLOBAL#MOAD_CEN_LAT=-8.4095078
>>>>>>>    NC_GLOBAL#STAND_LON=115.02
>>>>>>>   NC_GLOBAL#TRUELAT1=-8.4095001
>>>>>>>
>>>>>>> but obviously this is not being interpreted or used.
>>>>>>>
>>>>>>> Trying to use gdalwarp:
>>>>>>>
>>>>>>> lees-mbp:full Lee$ gdalwarp NETCDF:test_full.nc:HGT
>>>>>>> test_full_HGT.tif
>>>>>>> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
>>>>>>> ERROR 1: Unable to compute a transformation between pixel/line
>>>>>>> and georeferenced coordinates for NETCDF:test_full.nc:HGT.
>>>>>>> There is no affine transformation and no GCPs.
>>>>>>>
>>>>>>> Can anyone explain what I need to do?
>>>>>>>
>>>>>>> Thanks,
>>>>>>> Lee
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> gdal-dev mailing list
>>>>>>> gdal-dev at lists.osgeo.org
>>>>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> gdal-dev mailing list
>>>>>> gdal-dev at lists.osgeo.org
>>>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>
>>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131126/b3d455ba/attachment-0001.html>

More information about the gdal-dev mailing list