1

I downloaded MODIS product MOD13A1 for India. Then mosaiced the granules in GeoTIFF format using rasterio Python library. I have the appropriate shapefile to mask the mosaiced raster.

After mosaicing granules:

Mosaiced Raster

mosaic shape: (1, 12001, 14401)

crs: CRS({'proj': 'utm', 'zone': 35, 'ellps': 'GRS80', 'units': 'm', 'no_defs': True}) (found using mosaic.crs of rasterio library)

meta of mosaic file:

enter image description here

Shapefile I am using (converted to png to upload here):

enter image description here

CRS of shapefile: [epsg:4326]

I want to mask the mosaic using shapefile above. To bring them in same CRS I converted mosaic to EPSG:4326

Code to reproject mosaic:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'EPSG:4326'

with rasterio.open('mark_2.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

with rasterio.open('reproject_mark_2.tif', 'w', **kwargs) as dst:
    for i in range(1, src.count + 1):
        reproject(
            source=rasterio.band(src, i),
            destination=rasterio.band(dst, i),
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=transform,
            dst_crs=dst_crs,
            resampling=Resampling.nearest)

Code to mask mosaic using shapefile:

import fiona
import rasterio
import rasterio.mask

with fiona.open("india_shape_file/IND_adm0.shp", "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]

with rasterio.open('mosaic.tif') as src:
out_image, out_transform = rasterio.mask.mask(src, features,
                                                    crop=True)
out_meta = src.meta.copy()

But something is going wrong.

Reprojected mosaic:

enter image description here

Meta of reprojected file:

enter image description here

After masking:

enter image description here

Mosaicing Code:

import rasterio
from rasterio.merge import merge
from rasterio.plot import show
import glob
import os

dirpath = '37758/*'
search_criteria = 'MOD13A1_A2010001*.tif'
q = os.path.join(dirpath, search_criteria)

files = glob.glob(q)

mosaic_files = []

for fp in files:
    src = rasterio.open(fp)
    mosaic_files.append(src)

mosaic, out_trans = merge(mosaic_files)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
                 "height": mosaic.shape[1],
                 "width": mosaic.shape[2],
                 "transform": out_trans,
                 "crs": "+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs 
"})

with rasterio.open('mark_2.tif', "w", **out_meta) as dest:
    dest.write(mosaic)

Reprojection is wrong. Can anyone tell me what am I doing wrong ? I have tried the same in QGIS but result is same. Finally, I have to do this programmatically.

8
  • I doubt that your mosaic is really in UTM 35N, which is not a suitable CRS for India anyway. Could it be sinusoidal instead?
    – AndreJ
    Commented Nov 17, 2018 at 18:49
  • Individual granules were sinusoidal projection but after mosaicing it shows UTM. Here is the meta of single granule: {'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 2400, 'height': 2400, 'count': 1, 'crs': CRS({'proj': 'sinu', 'lon_0': 0, 'x_0': 0, 'y_0': 0, 'a': 6371007.181, 'b': 6371007.181, 'units': 'm', 'no_defs': True}), 'transform': Affine(463.31271652750013, 0.0, 4447802.078667, 0.0, -463.31271652750013, 4447802.078667)}
    – Sarvam
    Commented Nov 17, 2018 at 19:04
  • The corner coordinates of the affine transformation let me think that the mosaicing result is still in the same sinusoidal projection. The images are still 2400px squares. They would be distorted in any UTM projection.
    – AndreJ
    Commented Nov 17, 2018 at 19:10
  • @AndreJ I don't understand. If it is still in sinusoidal projection why is it showing otherwise. Also, after mosaicing why would image be 2400 px (or are you referring to granules ?). Could you please explain what is going wrong and what can be done ?
    – Sarvam
    Commented Nov 17, 2018 at 19:16
  • 1
    So you just wrote the UTM 35 in the metadata, but that did not reproject the raster to that. You better insert the sinu projection there.
    – AndreJ
    Commented Nov 17, 2018 at 19:47

1 Answer 1

1

In mosaicing code snippet out_meta has been updated manually.

out_meta.update({"driver": "GTiff",
             "height": mosaic.shape[1],
             "width": mosaic.shape[2],
             "transform": out_trans,
             "crs": "+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs 
"})

"crs": "+proj=utm" is setting the projection of mosaiced file to UTM when it should be sinusoidal.

Correct code will be:

"crs": "+proj=sinu +zone=35 +ellps=GRS80 +units=m +no_defs"
1
  • +zone=35 does not make sense for a sinusoidal projection. You should replace the ellipsoid with the a and b parameters you reported from a single MODIS image in your comment above. Or try to take the whole crs information from it.
    – AndreJ
    Commented Nov 18, 2018 at 7:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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