[R-sig-Geo] converting the projection from Lambert to WGS84 or NAD83

Wed Sep 20 21:50:38 CEST 2017

I am trying to reproject the eMODIS raster

  *   Projected Coordinate System:Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area
  *   Projection: Lambert_Azimuthal_Equal_Area
  *   Datum: D_Sphere_ARC_INFO

to the WGS84 or NAD83 18N

So that I can clip eMODIS raster files (more than 100) using a watershed boundary.
I attached the script here, but I still can not properly reproject the eMOIDS raster files.

# reading emodis data

# east US map tree phenology map
tmp<-raster("../../obs/eMODIS_East_SOSTAvg0114_v1/sostavg0114")

proj4string(tmp)
[1] "+proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

print(tmp)

> print(tmp)
class       : RasterLayer
dimensions  : 11556, 9174, 106014744  (nrow, ncol, ncell)
resolution  : 250, 250  (x, y)
extent      : 243000, 2536500, -2136500, 752500  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /work.taichi/Projects/Biscuit/rhessys/obs/eMODIS_East_SOSTAvg0114_v1/sostavg0114
names       : sostavg0114
values      : -115, 1000  (min, max)

# small watershed boundary map in NY state
bb<-raster("../../obs/gis/bb_soil_def_map.tif")
proj4string(bb)
[1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

print(bb)

> print(bb)
class       : RasterLayer
dimensions  : 598, 515, 307970  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 540217, 545367, 4648219, 4654199  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /work.taichi/Projects/Biscuit/rhessys/obs/gis/bb_soil_def_map.tif
names       : bb_soil_def_map
values      : 0, 5712  (min, max)

ex <- projectExtent(bb,  projection(tmp))

print(ex)

> print(ex)
class       : RasterLayer
dimensions  : 598, 515, 307970  (nrow, ncol, ncell)
resolution  : 30.12012, 14.94084  (x, y)
extent      : 6691756, 6707267, 9941473, 9950408  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

> ex
class       : RasterLayer
dimensions  : 598, 515, 307970  (nrow, ncol, ncell)
resolution  : 30.12012, 14.94084  (x, y)
extent      : 6691756, 6707267, 9941473, 9950408  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

c<-crop(tmp,  ex)

Error in .local(x, y, ...) : extents do not overlap

________________________________
From: Michael Sumner <mdsumner at gmail.com>
Sent: Tuesday, September 19, 2017 9:15 PM
To: Kyongho Son; RsigGeo
Subject: Re: [R-sig-Geo] converting the projection

Great thanks, that is good detail -but we still need the extents of both data to be able to see where the problem is.

You gave us the proj4string for both, but we also need the extent - you can print each object for a full summary eg.

print(bb)

print(tmp)

and with those anyone can reproduce enough (without any real data) to recreate your experience. (I am so used to working with this stuff I forget to remind others that print(theraster) is usually more than enough information to work through a reproducible example)

But, a guess is - unless the projection metadata is wrong, raster is probably telling you the truth - but it should be easy to see why there's a problem with the print output (essentially the proj4string and extent provide enough information for an experienced problem solver).

Cheers, Mike.

On Wed, 20 Sep 2017 at 11:05 Kyongho Son <kkyong77 at hotmail.com<mailto:kkyong77 at hotmail.com>> wrote:

Sorry and thanks, Mike
I can not share all codes and data since the original data are huge.
I copied the scripts that I used to reproject one of rasters and then make sure that they have same projection

 Then clip a large raster with a small raster. However, I got the error messages, the two raster maps are not overlap.

They should overlap if the projection was correct.

# east US map tree phenology map
tmp<-raster("../../obs/eMODIS_East_SOSTAvg0114_v1/sostavg0114")

proj4string(tmp)
[1] "+proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

# small watershed boundary map in NY state
bb<-raster("../../obs/gis/bb_soil_def_map.tif")
proj4string(bb)
[1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

ex <- projectExtent(bb,  projection(tmp))

> ex
class       : RasterLayer
dimensions  : 598, 515, 307970  (nrow, ncol, ncell)
resolution  : 30.12012, 14.94084  (x, y)
extent      : 6691756, 6707267, 9941473, 9950408  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

c<-crop(tmp,  ex)

Error in .local(x, y, ...) : extents do not overlap

*Kyongho Son*
Postdoctoral Fellow
Research Foundation of The City University of New York

________________________________
From: Michael Sumner <mdsumner at gmail.com<mailto:mdsumner at gmail.com>>
Sent: Tuesday, September 19, 2017 8:32 PM
To: Kyongho Son; RsigGeo

Subject: Re: [R-sig-Geo] converting the projection
Keep on list please.

Share code to demonstrate so we can see the problem.  Use projectRaster(a, b)  for alignment,  try crop(b,  ex,  snap = "out")  for inclusive crop.

And apologies I wrote clip instead of crop earlier.

Cheers,  Mike

On Wed, 20 Sep 2017, 07:42 Kyongho Son <kkyong77 at hotmail.com<mailto:kkyong77 at hotmail.com>> wrote:

I did the exact same way, but the output of the reprojection of 'a' raster did not properly align with 'b' raster.

*Kyongho Son*
Postdoctoral Fellow
Research Foundation of The City University of New York

________________________________
From: Michael Sumner <mdsumner at gmail.com<mailto:mdsumner at gmail.com>>
Sent: Tuesday, September 19, 2017 5:03 PM
To: Kyongho Son; r-sig-geo
Subject: Re: [R-sig-Geo] converting the projection

Try projecting just the extent
library(raster)
ex <- projectExtent(a,  projection(b))
clip(b,  ex)

where a,  b are your two raster.

Cheers,  Mike

On Wed, 20 Sep 2017, 06:18 Kyongho Son <kkyong77 at hotmail.com<mailto:kkyong77 at hotmail.com>> wrote:

Hi I have two raster data sets but they have two datum and projection.

I would like to convert one to the other since I would like to clip the one raster using the other raster.

they should have same projections.

proj4string(bb)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(SOSTavg)

 proj4string(SOSTavg)
[1] "+proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

Thanks,

Kyongho

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:R-sig-Geo at r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
<https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia+%3Chttps://maps.google.com/?q%3D203%2BChannel%2BHighway%2B%250D%2BKingston%2BTasmania%2B7050%2BAustralia%26entry%3Dgmail%26source%3Dg%3E&entry=gmail&source=g><https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway<https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
Kingston Tasmania 7050 Australia<https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>

--
<https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway<https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia+%3Chttps://maps.google.com/?q%3D203%2BChannel%2BHighway%2B%250D%2BKingston%2BTasmania%2B7050%2BAustralia%26entry%3Dgmail%26source%3Dg%3E&entry=gmail&source=g>
Kingston Tasmania 7050 Australia<https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia+%3Chttps://maps.google.com/?q%3D203%2BChannel%2BHighway%2B%250D%2BKingston%2BTasmania%2B7050%2BAustralia%26entry%3Dgmail%26source%3Dg%3E&entry=gmail&source=g>

--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]