7

I would like to create maps of a "tilted Earth", such as the following example where Australia is at the North Pole:

enter image description here http://i1048.photobucket.com/albums/s379/laskaris_mods/tilted_earth.jpg

I can reproject maps in image formats like .jpeg like this with the software I am using, Geocart 3. However, what I need to do is reproject elevation datasets in netcdf format (the ETOPO .grd data) in just this way, so that I have exact elevation data for my "tilted Earth".

Is there a GIS software that can do this?

Laskaris
  • 71
  • 2
  • Also asked on Cartotalk: http://www.cartotalk.com/index.php?showtopic=9050 where daan strebe reported that Geocart 3 does support TIFFs (you have to identify the TIFFs CRS in the software). – mkennedy Aug 09 '13 at 18:29
  • 1
    But I need to reproject the actual elevation data of the TIFF, and I don't think Geocart does that. I would also prefer to be able to do it with netcdf format. – Laskaris Aug 09 '13 at 18:37
  • Because the cell sizes/locations have changed, the cell values have to be resampled/recalculated but that should be happening. Why do you think it isn't? Reprojecting a point doesn't affect its elevation/z/H value. – mkennedy Aug 09 '13 at 18:49
  • You are right, it works for TIFF with Geocart. I was mistaken there. Thank you very much! Now, if I could do the same with netcdf, it would be even better. – Laskaris Aug 09 '13 at 19:02

1 Answers1

1

It can be done using python netCDF4, Projection4 library pyproj and numpy. Provided knowing the netcdf file CRS and tilted earth CRS, if EPSG codes are available it is super easy. The steps are

  1. First import netcdf file into python using netCDF4 library.
  2. Query the lat, long of netcdf file and store it as array.
  3. Define the current netcdf file CRS and required tilted earth CRS, if no CRS, the projection can be defined, for an example to convert between lcc custom projection into EPSG:4326, below code works. import pyproj
    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) lat2,lon2=[],[] for k, l in zip(lon,lat): lat1,lon1=pyproj.transform(proj1,proj2,k,l) lat2.append(lon1) lon2.append(lat1)
  4. The arrays lat2 and lon2 contains reprojected lat long values, use that values to recreate Netcdf file or tiff using python ogr library.

This answer may give a starter for steps 1 to 2. Or see this extended introduction for projection conversion using python library pyproj.

Community
  • 1
nish
  • 159
  • 7