Have LCC center terms, need PROJ.4 terms

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

Have LCC center terms, need PROJ.4 terms

hamish-2
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!

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

hamish-2
Hamish wrote:

> 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.


sorry, the attachment seems to have been stripped off by the ML. here 'tis:


 http://muck.otago.ac.nz/~hamish/sbss/lcc_grid_mismatch.png


Hamish

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

Andre Joost
In reply to this post by hamish-2
Am 09.02.2014 21:59, schrieb 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 ;

I guess that should be a LCC 1SP projection with these parameters:
+proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=0 +y_0=0
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs


> This trial and error gets me to within about 5% (see attatched)

 From the image it is hard to tell what grid line should be right or
wrong. For lon_0=-84, the 84W should be exactly vertical.

Can you share sample coordinates?

Greetings,
André Joost




_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

support.mn
In reply to this post by hamish-2
Hello,

after looking the picture you supplied it seems like
there is a rotation and scaling needed or also using
a different center setup for the projection might fix
that rotation:

1) try some other center lon or lat candidates

2) apply some rotation and scaling either using
a 7 parameter datum shift or

3) rotate and scale the resulting projection

It is most likely that any of those 3 corrections
will fix it but if the problem is due to wrong parameters
then number 1 would be the most desirable
solution.

Regards: Janne.

----------------------------------------------


Hamish [[hidden email]] kirjoitti:

> Hamish wrote:
>
> > 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.
>
>
> sorry, the attachment seems to have been stripped off by the ML. here 'tis:
>
>
>  http://muck.otago.ac.nz/~hamish/sbss/lcc_grid_mismatch.png
>
>
> Hamish
>


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

hamish-2
schrieb Hamish:

> > 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 ;

André wrote:
> I guess that should be a LCC 1SP projection with these parameters:
> +proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=0 +y_0=0
> +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Thanks, unfortunately I tried that & it's nowhere close.


> > This trial and error gets me to within about 5% (see attached)
[ http://muck.otago.ac.nz/~hamish/sbss/lcc_grid_mismatch.png ]

> From the image it is hard to tell what grid line should be right or
> wrong.

Sorry, I should have explained that: the black lines are the native
grid (contour lines created from the lat/lon raster grid values), and
the grey is an overlay of what my best-guess projection terms produce.

With a bit more trial and error using 1SP I've now got it matching a
little better, but it's not really satisfying to have to fudge it with
an random guess.

revised trial and error:
  http://muck.otago.ac.nz/~hamish/sbss/lcc_grid_mismatch2.png

-PROJ_INFO-------------------------------------------------
name       : Lambert Conformal Conic
proj       : lcc
a          : 6371229.000
es         : 0.0
f          : 0.0
lat_1      : 40
lat_0      : 33.05
lon_0      : -84.0
x_0        : 235000
y_0        : 0.000
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1.0

GRASS> g.region -gcl | grep center
center_long=-77.00767527
center_lat=39.85944826

but I'd really like to find the mathematical/code relationship between
the setup values I was given and what comes out the other end. Based
on other output grids I'm pretty happy that the lat/lon values in the
grids are ok, and it's the provided lcc terms which are the problem.


> For lon_0=-84, the 84W should be exactly vertical.

With +lon_0=-84 I need to set +x_0=235000 (or thereabouts) in order to
get the 84W meridian to be vertical; with +x_0 it is vertical at the
very left of the grid which according to the longitude grid values is
-87.018799. (the grid is 138x138 so I'm wildly guessing cell-center
registered not grid-conflence registered)

.. but +lon_0=-84 && +x_0=0 doesn't do it. I'd be surprised if the code
used a false easting, so I suspect the +lon_0 value really isn't 84W.
(and it's just on a sphere --nothing fancy, so I'd also be surprised
if there was any rotation or scaling going on Janne)


> Can you share sample coordinates?

sure, but I only have the grid and the lat/lon equivalents of the grid
and the angle-to-true-north grid data to go by, I don't have known-good
projected coordinates to test against beyond the center-lat/lon which
I can calculate. The x,y resolution "should" be 12000m, so the north,
east bounds in the following ascii grids are simply based on that
times the number of rows,cols (138x138):

  http://muck.otago.ac.nz/~hamish/sbss/xlat.txt.bz2
  http://muck.otago.ac.nz/~hamish/sbss/xlong.txt.bz2



thanks,
Hamish

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

Andre Joost
Am 12.02.2014 00:33, schrieb Hamish:

>
>    http://muck.otago.ac.nz/~hamish/sbss/xlat.txt.bz2
>    http://muck.otago.ac.nz/~hamish/sbss/xlong.txt.bz2
>



I checked it the other way round:

Changed the header of your files to:
> nrows        138
> ncols        138
> xllcorner    0
> yllcorner    0
> cellsize     12000
> NODATA_value -32768

which makes them readable for the GDAL AAIGrid importer.
Then created contour lines from the data, and overlayed it with a WGS84
grid.

The custom projection I came best with was

+proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=233925
+y_0=722000 +a=6371229 +b=6371229 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

the x_0 value was easy to get, just the distance between left margin and
the 84° meridian. For y_0, I measured the difference between WGS84 grid
and the countour lines grid. So, empirical too, but I think you can not
get much better from a 138x138 cell grid at that resolution.

Overlaying in QGIS I got this picture:

<http://home.arcor.de/andre.joost/WRF-LCC.png>

The cell values around the map center read (from line 75 and 76):

lon:  -77.068146/ -77.057312 / -76.927246 /-76.916199
lat: 40.32201 / 40.313641 / 40.214413 / 40.206055

So CEN_LAT and CEN_LON are just the average of the four middle cells,
and have nothing to do with the following projection parameters.

Just some more references on the WRF projection parameters:

http://mailman.ucar.edu/pipermail/wrf-users/2007/000698.html

http://www.mmm.ucar.edu/wrf/users/docs/user_guide/users_guide_chap3.html

HTH,
André Joost

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

Andre Joost
Am 12.02.2014 21:55, schrieb Andre Joost:

>
> The custom projection I came best with was
>
> +proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=233925
> +y_0=722000 +a=6371229 +b=6371229 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
>
> the x_0 value was easy to get, just the distance between left margin and
> the 84° meridian. For y_0, I measured the difference between WGS84 grid
> and the countour lines grid. So, empirical too, but I think you can not
> get much better from a 138x138 cell grid at that resolution.
>

Ok, now the analytical version ;-)

Put the center coordinates in a text file WRF4326.txt, and perform

cs2cs +init=epsg:4326 +to +proj=lcc +lat_1=39.338 +lat_0=39.338
+lon_0=-84 +k_0=1 +x_0=0 +y_0=0 +a=6371229 +b=6371229
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs WRF4326.txt >out.txt

on it. Subtract the result

595747.57 104990.12

from 828000 (that is 138x1200/2, the center of the raster grid)

and you get

x_0=232252. and y_0=723010.

Almost the same as I got empirical, taking the raster cell size of
12000m into acount.

HTH,
André Joost



_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Have LCC center terms, need PROJ.4 terms

Andre Joost
Am 12.02.2014 22:41, schrieb Andre Joost:

> Subtract the result
>
> 595747.57 104990.12
>
> from 828000 (that is 138x1200/2, the center of the raster grid)
>
> and you get
>
> x_0=232252. and y_0=723010.
>

Just to get it right:
the inversed values should go into the lower left corner values of the
ASCII Grid header:

nrows        138
ncols        138
xllcorner    -232252.
yllcorner    -723010.
cellsize     12000
NODATA_value -32768

and with that no false Easting or Northing is necessary for the projection:

+proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=0 +y_0=0
+a=6371229 +b=6371229 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Greetings,
André Joost

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj