3

I have two numpy meshgrids (X, Y both 2D) in Lambert Conformal conic with the following information:

enter image description here

I'd like to directly (or indirectly if not possible) project them to UTM. I have both pyproj and gdal installed, but wasn't able to figure out how to use them. I found solutions for converting from LCC to WGS84 here and here, but couldn't use either of them since I don't have all the information needed.

Even if somehow I can convert from LCC to WGS84, then I can use Python UTM package and utm.from_latlon to get UTM out of it.

I found the following line of code from here:

awips221 = Proj(proj='lcc', R=6371200, lat_1=50, lat_2=50,lon_0=-107, ellps='clrk66')

And then I can use transform from pyproj. In that case, is lat_1 and lat_2 both 25.0 in my case? What would be lon_0 and ellps?

Using gdalinfo we can see lCC and WGS84 projections for the center and four corners of the grid in the picture below:

Projection description using <code>gdalinfo</code>

Improve this question
13
  • 1
    Try lon_0 -95 and ellps wgs84.
    – AndreJ
    Aug 19 '15 at 4:18
  • Following your suggestion, I use the transform function (from pyproj import transform). Then dine lccProj using lccProj = Proj(proj='lcc', R=6371229, lat_1=25, lat_2=25,lon_0=-95, ellps='wgs84') and utmProj using utmProj = Proj(proj='utm', zone=10, ellps='WGS84'). Finally, I transform them with XX_UTM, YY_UTM = transform(lccProj, utmProj, XX, YY). However, I am not getting meaningful results. I am using zone=10 for no good reason here! Can that be the cause or I am missing something else?
    – ahoosh
    Aug 19 '15 at 17:14
  • From this link, using lccProj defined the the previous comment, I can also do lon, lat = lccProj(XX, YY, inverse=True) and get longitude and latitude. But these value are not correct for some reason.
    – ahoosh
    Aug 19 '15 at 17:41
  • 1
    Using GDAL cs2cs, I get -133.45901182°E 12.26975675°N for lcc bottom left. I'm not sure whom to trust.
    – AndreJ
    Aug 19 '15 at 18:59
  • 1
    Maybe. I am using this projection string: +proj=lcc +lat_1=25 +lat_0=25 +lon_0=-95 +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. According to stackoverflow.com/questions/26452972/… you should be able to use it with pyproj.Proj
    – AndreJ
    Aug 19 '15 at 19:17
2

I'm not sure what NOAA thinks are the right coordinates, but I have no problem loading the file into QGIS, or reprojecting it to WGS84 with gdalwarp:

enter image description here

QGIS uses this custom projection string:

+proj=lcc +lat_1=25 +lat_2=25 +lat_0=25 +lon_0=265 +x_0=0 +y_0=0 +a=6371229 +b=6371229 +units=m +no_defs

You can use the same string with pyproj.Proj().

where you can safely replace lon_0=265 with lon_0=-95.

Note that the extent in degrees is not the same as reprojecting the corner coordinates, because the picture is heavily bended.

I'm not sure why you want to transform the data to UTM. It crosses several UTM zones (7N to 21N). You can take the middle one (15N), and experience distortions on the sides, or cut the picture along the UTM zone borders, and make calculations within each UTM zone separately.

Improve this answer
2
  • 1
    I was going to say that this projection is very strange. The central meridian crosses the +180 threshold and the spheroid parameters are unusual. It looks close to the GRS80 authalic sphere, but it's not quite there. Nice job.
    – Mintx
    Aug 20 '15 at 19:28
  • This works great with good accuracy. For example, lower left corner which is [-4228464.497,-835237.761] in 'lcc' should result in [-133.474506,12.162367] based on photo from gdalinfo in my question. If I use what you suggested with pyproj.Proj(), I get [-133.472937, 12.162813] which is just 177 m apart using Geod function. I will do my calculation within each UTM zone separately for sure. Thanks a bunch.
    – ahoosh
    Aug 20 '15 at 19:42
0

I tried this string in pyproj:

 Proj(proj='lcc', R=6371200, lat_1=25, lat_2=25,lon_0=-95, ellps='clrk66')

and my results were strange. when I changed the lat field names:

 Proj(proj='lcc', R=6371200, lat_0=25, lat_1=25,lon_0=-95, ellps='clrk66')

I got the results I was expecting.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.