<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Use pyresample to make a projected image2&nbsp;&nbsp;Read the lons/lats from the MYD03 file3&nbsp;&nbsp;get the map projection from corners.json4&nbsp;&nbsp;Use pyresample to define a new grid in this projection5&nbsp;&nbsp;resample the longitudes on this grid6&nbsp;&nbsp;replace missing values with floating point nan7&nbsp;&nbsp;Plot the image using cartopy  Use pyresample to make a projected imageIn the cartopy_mapping_pyproj notebook we stored projectioncoords in a json file called corners.json.  This notebookreads that information back in to plot lats/lons on a map<jupyter_code>import a301
import json
from a301.utils.data_read import download
import a301
import pprint
import shutil
from pyhdf.SD import SD, SDC
import json
import pprint
import cartopy
from pyresample import kd_tree

read_data=False
if read_data:
    filename_M3='MYD03.A2013222.2105.006.2013223155808.hdf'
    download(filename_M3)
    for filename in [filename_M3,filename_M2]:
        local_file = Path.cwd() / Path(filename)
        to_file = a301.data_dir / Path(filename)
        print(f'copy {local_file} to {to_file}')
        shutil.copy(local_file,to_file)
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy
from pathlib import Path
import pprint
import numpy as np
import pdb
#<jupyter_output><empty_output><jupyter_text>Read the lons/lats from the MYD03 file**substitute your filename**<jupyter_code># Read the lats and lons from the MYD03 file
filename_M3='MYD03.A2013222.2105.006.2013223155808.hdf'
m3_path= a301.data_dir / Path(filename_M3)
print(f'reading {m3_path}')
m3_file = SD(str(m3_path), SDC.READ)
lats = m3_file.select('Latitude').get()
lons = m3_file.select('Longitude').get()<jupyter_output>reading /Users/phil/repos/a301_code/data/MYD03.A2013222.2105.006.2013223155808.hdf<jupyter_text>get the map projection from corners.jsonGet the map  projection and extent from corners.json<jupyter_code>json_file = a301.data_dir / Path('corners.json')
with open(json_file,'r') as f:
    map_dict=json.load(f)
pprint.pprint(map_dict)<jupyter_output>{'extent': [-1285873.5967137816,
            1561347.9917805532,
            -1179100.5032042824,
            1297248.5261361937],
 'lats': [32.13645206898284,
          28.687374622563773,
          45.73346985640787,
          50.510827489422674,
          32.13645206898284],
 'lons': [-104.77089390290801,
          -129.005397891393,
          -138.038848796623,
          -107.001718605882,
          -104.77089390290801],
 'proj4_params': {'datum': 'WGS84',
                  'ellps': 'WGS84',
                  'lat_0': 39.59910106367865,
                  'lon_0': -121.4048713497655,
                  'proj': 'laea',
                  'x_0': 0.0,
                  'y_0': 0.0},
 'proj4_string': '+datum=WGS84 +ellps=WGS84 +proj=laea '
                 '+lon_0=-121.4048713497655 +lat_0=39.59910106367865 +x_0=0.0 '
                 '+y_0=0.0 +no_defs',
 'xcoords': [1561347.9917805532,
             -744961.1366254934,
             -1285873.5967137816,
             1019738.9399581843,
   [...]<jupyter_text>Use pyresample to define a new grid in this projection<jupyter_code>from pyresample import load_area, save_quicklook, SwathDefinition
proj_params = map_dict['proj4_params']
swath_def = SwathDefinition(lons, lats)
area_def=swath_def.compute_optimal_bb_area(proj_dict=proj_params)
dir(area_def)<jupyter_output><empty_output><jupyter_text>resample the longitudes on this grid<jupyter_code>fill_value=-9999.
area_name = 'modis swath 5min granule'
image_lons = kd_tree.resample_nearest(swath_def, lons.ravel(),
                                  area_def, radius_of_influence=5000, 
                                      nprocs=2,fill_value=fill_value)
print(f'\ndump area definition:\n{area_def}\n')
print((f'\nx and y pixel dimensions in meters:'
       f'\n{area_def.pixel_size_x}\n{area_def.pixel_size_y}\n'))<jupyter_output>dump area definition:
Area ID: laea_otf
Description: On-the-fly laea area
Projection: {'datum': 'WGS84', 'ellps': 'WGS84', 'lat_0': '39.59910106367865', 'lon_0': '-121.4048713497655', 'proj': 'laea', 'x_0': '0.0', 'y_0': '0.0'}
Number of columns: 1489
Number of rows: 2244
Area extent: (-1269560.846045296, -1177850.7780355075, 1532203.80246403, 1265820.4615973115)


x and y pixel dimensions in meters:
1881.641805580474
1088.980053312308<jupyter_text>replace missing values with floating point nan<jupyter_code>nan_value = np.array([np.nan],dtype=np.float32)[0]
image_lons[image_lons< -9000]=nan_value<jupyter_output><empty_output><jupyter_text>Plot the image using cartopy<jupyter_code>crs = area_def.to_cartopy_crs()
ax = plt.axes(projection=crs)
ax.coastlines()
ax.set_global()
plt.imshow(image_lons, transform=crs, extent=crs.bounds, origin='upper')
plt.colorbar();
crs.globe.to_proj4_params()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup