[gdal-dev] WRF NetCDF import

Lee Eddington lee.w.eddington at gmail.com
Tue Nov 26 02:32:13 PST 2013 

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.

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
>>>>>
>>>>
>>>>
>>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131126/eaf1c354/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: read_wrf_netcdf_U10.py
Type: text/x-python-script
Size: 1837 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131126/eaf1c354/attachment-0001.bin>

More information about the gdal-dev mailing list