[R-sig-Geo] converting the projection from Lambert to WGS84 or NAD83
Francisco Rodriguez Sanchez
f.rodriguez.sanc at gmail.com
Thu Sep 21 12:58:41 CEST 2017
Hi Kyongho,
There seems to be something strange with your eMODIS raster, which I
understood represents something in the eastern US. There might be
something wrong with the coordinates or the original projection
information. When I reprojected it to global lonlat to see which area it
covers, it actually fell in Antartica! (see
https://i.imgur.com/wjTBt8E.png). So, I'm not sure if I got it
correctly, but I think that's not what you expect. That would explain
the error that the extents do not overlap.
Reproducible example follows:
library(raster)
## Create raster with same extent and projection
emodis <- raster(xmn = 243000, xmx = 2536500, ymn = -2136500, ymx = 752500,
resolution = 10000,
crs = "+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",
vals = 0)
emodis
#> class : RasterLayer
#> dimensions : 289, 229, 66181 (nrow, ncol, ncell)
#> resolution : 10000, 10000 (x, y)
#> extent : 243000, 2533000, -2137500, 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 : in memory
#> names : layer
#> values : 0, 0 (min, max)
## Project to lonlat to see where it falls
emodis.lonlat <- projectRaster(emodis, crs = "+proj=longlat +datum=WGS84")
library(maptools)
data(wrld_simpl)
plot(wrld_simpl)
plot(emodis.lonlat, add = TRUE)
# See map at 
Hope this helps!
Cheers,
Paco
El 20/09/2017 a las 21:50, Kyongho Son escribió:
> 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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Dr Francisco Rodriguez-Sanchez
Integrative Ecology Group
Estacion Biologica de Doñana (CSIC)
Avda. Americo Vespucio 26
E-41092 Sevilla (Spain)
http://bit.ly/frod_san
More information about the R-sig-Geo
mailing list