# Conversion of Latitude/Longitude to different map projection

I have two m*n matrices coming from a WRF-ARW simulation: one exclusively with longitude values and one with only latitude values. For both matrices I know the originatin projection (lambert conic conformal, with true latitudes and parameter set by me).

I would like to convert them to a regular plate carree projection (the one employed by Google Earth, to be precise): is it possible to do so? If yes, which would be the best instrument and how, exactly?

I'm currently having at hand R (where I've and may have both QGis and MATLAB (but without any toolbox)

Thanks for any help!

You can add the matrices as ASCII Grid file to QGIS, and create contours from the lat and lon values. This will create a degree grid like this:

In this mailing list topic I described how to set up a suitable projection for it:

http://osgeo-org.1560.x6.nabble.com/Have-LCC-center-terms-need-PROJ-4-terms-td5102799.html

EDIT

Sample data can be found in the mailing list thread. The known projection parameters are:

``````TRULAT1 = 39.338
TRUELAT2 = 39.338
STAND_LON = -84.
``````

This gives the following proj definition string:

``````+proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
``````

No ellipsoid is given, so we take WGS84 as default.

The extent information inside the data reads:

``````north: 1656000
south: 0
east: 1656000
west: 0
rows: 138
cols: 138
``````

which can be transformed to the GDAL-known AAIGrid header:

``````nrows        138
ncols        138
xllcorner    0
yllcorner    0
cellsize     12000
NODATA_value -32768
``````

To get the x_0 and y_0 values, we take the lower left values from the data (last line,first value): `-86.433044 33.062443` into a file `ll.txt` and convert those to lcc using cs2cs:

``````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 +datum=WGS84 +units=m +no_defs ll.txt >>lcc84.txt
``````

Resulting to `-228512.96 -694637.98 0.00`.

Note that the cells are 12km wide, and the cell value is for mid cell. So we have to invert the values and add 6000m to get the real lower left values:

``````234513.  700638.
``````

so the complete proj string is:

``````+proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=234513 +y_0=700638 +datum=wgs84 +units=m +no_defs
``````

Another approach is to look up the lon_0 value in the xlong table. It is in every line at column 20. So the x_0 is 19.5 * 12000 = 234000. Looking for lat_0 in the xlat table at the same column returns row 59 from the bottom, so the correct y_0 value is 58.5 * 12000 = 702000. The improved proj string is:

``````+proj=lcc +lat_1=39.338 +lat_0=39.338 +lon_0=-84 +k_0=1 +x_0=234000 +y_0=702000 +datum=wgs84 +units=m +no_defs
``````

A more visual method is to translate the asc files into xyz files, and load those into QGIS as delimited text:

``````gdal_translate -of XYZ xlong.asc xlong.xyz
gdal_translate -of XYZ xlat.asc xlat.xyz
``````

using the first lcc proj string with x_0=0 and y_0=0. Then look up the coordinates of lon_0 and lat_0, they are 234000 and 702000; which approves the last mentioned values for x_0 and y_0.

• Following your suggestion, I am unable to do the re projection using pyproj library. It is giving wrong lat long with units in centimeters. I defined the wrf arw projection in python as follows, `tc={'proj':'lcc','width':width_meters,'height':height_meters,'lat_0':cen_lat,'lon_0':cen_lon,'lat_1':truelat1,'lat_2':truelat2}; proj1=pyproj.Proj(tc); proj2_out='+init=EPSG:4326'; proj2=pyproj.Proj(proj2_out); lat1,lon1=pyproj.transform(proj1,proj2,x,y);` any idea?
– nish
Nov 25, 2014 at 5:55
• width and height are no parameters of the lcc projection. Maybe that spoils the transformation. x_0 and y_0 are missing in your definition. Nov 25, 2014 at 6:38
• Can you please provide formula to calculate x_0, y_0. I tried to follow your analytical version solution as per this but not giving the result as mentioned. What I did is instead of`cs2cs`, used pyproj, `proj1=pyproj.Proj("+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") proj1=pyproj.Proj(tc) proj2_out='+init=EPSG:4326' proj2=pyproj.Proj(proj2_out)` but giving coordinates of (597909.4356527944, 60082.23310022025)
– nish
Nov 25, 2014 at 9:08
• For which input values? Nov 25, 2014 at 9:14
• See my extended answer. Nov 25, 2014 at 20:36

Below is a work around method using python and bunch of libraries such as netCDF4, numpy and shapefile. With Anaconda, these library installations are very much easy. The steps of the method are

1. Import WRF ARW output into python, by python netCDF4 library.
2. Query WRF ARW output variable XLAT, XLONG into numpy array of latitude and longitude using python numpy library.
3. loop into the latitude and longitude, create point shape file using python shape file library.
4. Open the shape file in QGIS and specify the required projection such as EPSG 4326. Now the WRF ARW output is in projected shape file formate to be able to converted into any form such as kml, usable in Google earth.

Here is an example

1. To import WRF ARW output into Python. ``` import shapefile as shp from netCDF4 import Dataset import numpy as np wrfoutput = 'wrfout_d01_2014-06-05_00:00:00' fh = Dataset(wrfoutput, mode='r') ```
2. Query WRF ARW output variable XLAT, XLONG. ``` lons = fh.variables['XLONG'][:] lats = fh.variables['XLAT'][:] #arrays are in 3D, reshape it lonsX0=lons.reshape((105,96)) #into 2D lonsX=lonsX0.reshape((10080)) #into 1D latsX0=lats.reshape((105,96)) latsX=latsX0.reshape((10080)) IN=np.array(range(10800)) #an index array for point shape file ```
3. Now loop into the 1D coordinate arrays (latsX,lonsX) and create shape file. ``` w = shp.Writer(shp.POINT) w.field("fo","N",10) for j,k in enumerate(lonsX): w.point(k,latsX[j]) #write the geometry w.record(IN[j]) out_file='WRF_output.shp' ```
4. Open it in QGIS and specify the required CRS and you are done.