Difference between revisions of "Building Surface Topography Files"

#!/usr/bin/env python3.0
#-*- coding:Utf-8 -*-

from    numpy import *
import  numpy                 as        np
import  matplotlib.pyplot     as        mpl
import  math
import xarray as xr

# 1. WE GET THE SURFACE TOPOGRAPHY DATA
nc = xr.open_dataset('surface_mars.nc',decode_times=False) # can be any netcdf file (e.g. start/startfi.nc files)

# 2. WE READ THE VARIABLES
lat=nc['latitude']
lon=nc['longitude']
albedo=nc['albedo']
thermal_inertia=nc['thermal']
elevation=nc['zMOL']

# BELOW THE VARIABLES WE WANT TO UPDATE
new_elevation = np.empty((len(lat),len(lon)))
new_thermal_inertia = np.empty((len(lat),len(lon)))
new_albedo = np.empty((len(lat),len(lon)))

# LOOP TO MODIFY THE VARIABLES
for i in range(0,len(lat),1):
    for j in range(0,len(lon),1):
        new_elevation[i,j]=max(elevation[i,j],-3.9) # here you put whatever you want ; in this example, we take the (MOLA) topography of present-day Mars and modify it to fill the topographic depressions with an ocean at an altitude of -3.9km
        new_albedo[i,j]=albedo[i,j] # here you put whatever you want
        new_thermal_inertia[i,j]=thermal_inertia[i,j] # here you put whatever you want

# SANITY CHECK PLOTS

fig1 = mpl.figure(1)
mpl.contourf(lon,lat,elevation)
mpl.xlabel('longitude (deg)')
mpl.ylabel('latitude (deg)')
mpl.show()

fig2 = mpl.figure(2)
mpl.contourf(lon,lat,new_elevation)
mpl.xlabel('longitude (deg)')
mpl.ylabel('latitude (deg)')
mpl.show()

# CREATE THE NEW NETCDF TOPOGRAPHY FILE

nc['albedo'].values = new_albedo
nc['zMOL'].values = new_elevation
nc['thermal'].values = new_thermal_inertia
nc.to_netcdf('new_topography.nc')