Conversion of Latitude/Longitude to different map projection

Question

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!

Answer 1

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
MOAD_CEN_LAT = 39.33799
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.

Answer 2

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.