[Proj] Have LCC center terms, need PROJ.4 terms

[Hamish]

Hi,

I'm working with some weather model (WRF) output and I'll like to
import the results into proj4-aware software (gdal or grass gis).

The LCC model domain specification takes center of region lat/lon as
input, not +lat_0, +lon_0. The values I've been given for +lat_1 and
+lat_2 are equal, and near to +lon_0. That doesn't look quite
right but the model seems happy enough with it:


the values I was given:


CEN_LAT = 40.26405f ;
CEN_LON = -76.99222f ;
TRULAT1 = 39.338f ;
TRUELAT2 = 39.338f ;
MOAD_CEN_LAT = 39.33799f ;
STAND_LON = -84.f ;


What I do have of value is two of the raster grids put out by the
model which contain the lat and lon of each grid point, and these
seem correct. I also have two raster outputs for convergence angle,
ie the angle between grid north and geographic north at each grid
point. I know that +x_0 and +y_0 are 0, and that +ellps is a sphere,
the WRF code lists about 4 different radii for +a, but that's a
secondary issue. Grid cell size is down to 4km, so I'm not overly
concerned with datum issues on the order of 100m.


What I'd like to do is to take the center-of-region LCC definition
and solve for the +lon_0, +lat_0, +lat_1, +lat_2 terms so I can
process the model's output grids in their native projection.

I could treat the two lat,lon grids as source data to run a manual raster
reprojection and interpolate into a proj4-know projection, but that's
computationally expensive and lossy.

So far I've gotten close by choosing the bottom-left grid point of
the lat and lon grids, and reading off those values to use for the
+lat_0 and +lon_0 terms. Then by educated-guessing I picked some
+lat_1 and +lat_2 terms and iteratively have gotten pretty close to
being able to overly a geographic grid on top of contour lines from
the lat and lon raster grids. Of course the terms I found are no
where near those I was provided with or that I can see in the model's
domain setup files.


my guess:

-PROJ_INFO-------------------------------------------------
name       : Lambert Conformal Conic
proj       : lcc
a          : 6371229.000
es         : 0.0
f          : 0.0
lat_0      : 33.0624
lat_1      : 20
lat_2      : 45
lon_0      : -87.0188
x_0        : 0.0000
y_0        : 0.0000
-PROJ_UNITS-------------------
-----------------------------
unit       : meter
units      : meters
meters     : 1.0


This trial and error gets me to within about 5% (see attatched), but
I'd really like to do the job properly and calculate the exact terms.
I have access to the Fortran code used in the model to make the lat
and lon grids, could I use them to also derive the proj4 terms?


thanks for any ideas or advice on how to proceed.


Hamish

ps- sorry for any crazy formatting, yahoo mail's web app is broken yet again and won't let me change to ascii or have any control over the formatting. grumble!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20140209/507deac9/attachment.htm