On Tue, Nov 26, 2013 at 8:32 AM, Lee Eddington <lee.w.edding...@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 < > rveci...@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.edding...@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.edding...@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....@gmail.com> wrote: >>>> >>>>> Forwarding to the list and others... >>>>> >>>>> ---------- Forwarded message ---------- >>>>> From: Roger Veciana i Rovira <rveci...@gmail.com> >>>>> Date: Mon, Nov 25, 2013 at 9:35 AM >>>>> Subject: Re: [gdal-dev] WRF NetCDF import >>>>> To: Etienne Tourigny <etourigny....@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....@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.edding...@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@lists.osgeo.org >>>>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >>>>>>> >>>>>> >>>>>> >>>>>> _______________________________________________ >>>>>> gdal-dev mailing list >>>>>> gdal-dev@lists.osgeo.org >>>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >>>>>> >>>>> >>>>> >>>>> >>>> >>> >> > > _______________________________________________ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev >
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev