texts: Error using tmap on large sf object in R

texts:
I'm trying to use Resolve ecoregion data to map the biomes of the African continent. The shapefile is global, and there are over 800 polygons, as each polygon represents a unique ecological area (whereas there are only seven terrestrial biomes, biomes are the larger spatial unit and are composed of similar ecological areas). Perhaps because of the data size or number of polygons, I am having a difficult time producing a map using tmap.
I first cropped the shapefile to the African continent:
resolve = st_read("Data/Ecoregions2017/Ecoregions2017.shp")

Africa = st_read("Data/Africa_SHP/Africa/Africa.shp")
st_crs(Africa) = crs(resolve)

resolveAfrica =st_intersection(st_make_valid(resolve), st_make_valid(Africa))

But when I try to then map this cropped shapefile:
tm_shape(Africa) + tm_borders() +
  tm_shape(resolveAfrica) + tm_fill("BIOME_NAME")

R returns an error:
Error in vapply(lst, class, rep(NA_character_, 3)) :
  values must be length 3,
 but FUN(X[[56]]) result is length 2

Does anyone know what steps I can take to address this error?


texts: Selecting multi features on map based on Datagrid row/rows selection?

texts:
I was trying to select feature depending on the selection change  of the datagrid. Each row of the datagrid consists of the feature OID and feature Name etc.

While selecting I am seeing that more features are selected in the map than I selected in the row. While investigating the issue, I found that a feature OID is unique through out its own attribute Table but it may also be available in different feature Attribute Table.

Here is the code #

foreach (DataGridViewRow row in dataGrid.SelectedRows)
{       grdrows.Add(Convert.ToInt32(dataGrid.SelectedRows[iCount].Cells[0].Value));
 iCount++;
 }
            int[] ids = grdrows.ToArray();
            map.ClearSelection();
            if (dataGrid.SelectedRows.Count > 0)
            {
              foreach (DataGridViewRow row in dataGrid.SelectedRows)
               {
                   int objectID = Convert.ToInt32(row.Cells[0].Value.ToString());
                   string featureName = row.Cells[1].Value.ToString();

                    while (counter != map.LayerCount - 1)
                     {
                       IFeatureLayer featureLayer =
                      (IFeatureLayer)map.get_Layer(counter);

                       ISelectionSet pSelSet =
                       ((IFeatureSelection)featureLayer).SelectionSet;
                       pSelSet.AddList(ids.Count(), ref ids[0]);
                      counter++;
                    }
            }
        }

Then I have added this if but it did not work. The problem is that post adding the if statement multiple row selection does not select multiple features. It selects only the last feature. I need help.
I have put the selection code in this .

if ( featureLayer.FeatureClass.AliasName == featureName )


texts: QGIS 3.12 print composer - how to rotate map but keep grid unrotated?

texts:
I am creating a map using QGIS 3.12 and wish to rotate the map and add a grid using print composer, but I want to keep the grid unrotated.

Is this possible or does the grid only rotate with the map?


texts: Converting many .dxf files to one file geodatabase

texts:
I am using ArcView and ArcCatalog 9.2 with the Data Interoperability extension and I have 159 .dxf files that are comprised of 6 files; a dessin DAO, an annotation DAO, a multipatch DAO; a point DAO, a polygon DAO, and a polyline DAO.

It's all in French and I don't know the english words for some of the tools. So I apologize for any unclarity.

They create a parcel layer with some building foot prints along with other features. Using a tool (Eclater un DAO; translation: Explode a DAO) in ArcCatalog that converts DAO to GDB works, but I will have to do it 159 times getting 159 different GDB when all I want is for them to be together.

Will I need to combine the 159 .dxf files into one first? If so, how do I do that without CAD?


texts:
We are using QGIS in our workgroup, and have recently updated to 1.8.0 Lisboa.
To share our file structure we are using egnyte similar to dropbox.
We all use Mac OS, and have used QGIS for years, and egnyte for several months now without any problems.

After updating to 1.8.0 we have experienced a frightening problem: The shx-files from a vast amount of our shape-folders are disappearing mysteriously, leaving the data useless. Fortunately we have been able to recover most of the shx-files from back ups, but obviously we can not rely on that for the future. The shx-files that disappear seem to be from folders we have recently worked with, but I'm not able to see any other tendencies in what is gone and what is still intact.

Very grateful for any help on this issue. All the best


texts:
Is there any tool that can help me with my problem?

I have polygons looking like this :

And what I want, is to simplify or gain convex boundary of such polygons.

I tried Simplify Polygon in ArcMap, Simplify Geometryin QGIS and v.clean in Grass using: snap lines/boundaries, remove small angles at nodes.

I tried these tools also on polygons converted to lines.

Some of these tools worked, but my layer was changed too much.
I need something like inteligent convex boundary.

Do you have any ideas ?


texts:
I'm working with Sentinel-2 data in QGIS and after opening the product I zoom several times to have a look to the image and the program crashes. I'm using QGIS 2.18.10 on a Windows 10 mahcine. Any idea about this error?

M

Here some screenshots of the error:


texts: Kernel Density (KDE) Issues and GME Errors: "Failed to create raster dataset" and "Arithmetic operation resulted in an overflow."

texts:
I'm trying to create 50% and 95% core and home ranges in GME and am running into some problems.

First, when I try to to create my KDE raster, occasionally I get a "Failed to create raster dataset" error. Sometimes I get it and sometimes I don't. As you can see here, all I did was change the name and it worked the second time. Anyone know what might be going on? My settings for the KDE raster are: Bandwidth = PLUGIN and Cell size = 30.

Second, if I do successfully create a KDE raster and try to create my home range raster using the ISOPLETH tool, I always get an error saying that the "Arithmetic operation resulted in an overflow". For the home range raster I've used multiple variations for the quantile (i.e., 0.5, 0.95, c(0.5, 0.95), 50, 95) and none of them have worked. Anyone know what might be going on?

Alternatively, I've tried to do the core and home range in ArcGIS 10.3.1 using the Kernel Density Tool (Spatial Analyst Tools). To get the 50% and 95% range, I need to manually reclassify the layer, correct? I've reduced the classes from 9 to 2 and manually set the breaks at 50% and 95%. Unfortunately, it's not producing a smooth home range polygon and has been producing a square for the 95% range. What am I doing wrong? Also, where can I find the actual metric area? While the graphics are nice, it's really the size of the actual area that I'm most interested in; I'm going to be comparing home range size between genders.


texts: Opening a popup on Marker click when Markers are in array using leaflet maps

texts:
Hey guys I have a set of Markers that are getting created dynamically on a .each loops like so.

$.each(data, function(index, element) {
    markerArr[element.id] = L.marker(map.unproject([element.LocationX,element.LocationY],map.getMaxZoom()-4)).addTo(map);
});

My question is how to open the popup for each popup on click and show the unique ID inside the popup?

I am able to do this if I  were to declare each Marker with a unique name like so.

var marker1 = L.marker([158,395]).addTo(mapl);
var marker2 = L.marker([158,410]).addTo(mapl);
marker1.bindPopup("Kanye West");
marker1.on('click', function (e) {
    this.openPopup();
});
marker2.bindPopup("50 Cent");
marker2.on('click', function (e) {
    this.openPopup();
});

But I need to be able to open all Markers that are in the array instead.

Tried something like this with no luck

markerArr[this].on('click', function (e) {
    this.openPopup();
});


texts: Setting the correct grids in QGIS 3.8 map composer?

texts: My area lies between E 92.25-92.50 and N 24.50-24.75. When I tried to draw the coordinates on QGIS 3.8 map composer, it showed the co ordinates as E 8.00-E 24.000 and N 0.000-16.000. I used the correct coordinates and projections but the problem still persists. I took "Decimal with suffix" and used "cm" or "mm" as the interval units. But that should not be a problem because the layer and the project are in the same projections. Please show with illustration how to fix this problem. QGIS grid outines are attached.


texts: If/Then Statement with LIKE using the Field Calculator

texts:
I cannot get the below to work with a LIKE variable within the If/Then statement. Is there something that I am missing?

Prelogic Script

dim y

If [Categories] LIKE '%Orange%' Then
y = "Medium"
Else y = [Categories]
End If

Status =
y


texts: Are special interpolation techniques required for weather data?

texts:
I've downloaded a series of weather observations from the local weather agency - it comprises around 1,000 points spread non-uniformly across Australia, with a greater density of points located near populated areas.

I'm interested in displaying mean temperatures and monthly rainfall (separately) as smooth surfaces, so I used the Kriging function in Spatial Analyst (with the default values) to calculate this grid:

[edit]

The grids will be used for visualisation purposes only. In a web map I'll display the grids to give context, but will only allow identify on the actual weather stations, in which case I'll show the original values. I guess this means the accuracy of the grids isn't of paramount importance.

My questions:

Which of the parameters on the Kriging tool do I need to understand, in order to create meaningful interpolations of rainfall, temperature and other weather factors?
Do I need to account for effects of topography (eg, orographic rain effect) or does the kriging algorithm handle this?
is the Geostatistical Analyst a better tool for this type of analysis? (I believe it's designed to handle biases due to geography)
anything else I should be aware of to create a meaningful result?

Thanks


texts: Generating report or table for buffer zone using ArcGIS Desktop?

texts:
I have a polygon feature class.

I create a buffer for this feature with a distance being 2 km.

Now the buffer output intersect with features of the original polygon.

So say for instance from the orginal polygon id 1 it creates a buffer polygon with id 1a.

Now I want a report or table of all the features the 1a intersects with the original polygon.

Any suggestions?


texts:
I want to extract every row from a .shp from a raster (Extract by Mask).
With the following code i get 4 output raster, wich concludes every row.

import arcpy
arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = 1
# path name variables
clipFeat = r'../all.shp'
output = r'/output'
raster = r'/Karte1.jpg'
# Name of OID field
fieldName = "FID"
cursor = arcpy.SearchCursor(clipFeat)
for poly in cursor:
    mask = str(poly.FID)
    # Create feature layer of current mask polygon
    whereClause = "poly.FID = poly.FID "
    arcpy.MakeFeatureLayer_management(clipFeat, 'currentMask',
    whereClause)
# Save the clipped raster
outRaster = arcpy.sa.ExtractByMask(raster, 'currentMask')
outRaster.save(output + "\\" + mask + ".tif") # Create your
own output name
arcpy.Delete_management('currentMask')
if arcpy.Exists('currentMask'):
    arcpy.Delete_management('currentMask')
    del cursor

If I type "whereClause = "\"FID\" = 0", I Got four files with the export of row 0.

My goal is that i got 4 differnt output, one für FID=0, FID=1, FID=2 and so on.


texts: SQL Server update statement using .NET made function does not run

texts:
In SQL Server I have a table called Track containing metadata of each track and a spatial representation of it: the field Geo, which is a LineString of the geography type.
In Visual Studio I developed a function MakeLineString that constructs those spatial representations based on GPS points in a table. The input for the function is the field TrackId of the table Track. I imported the function in SQL Server and it runs well in itself, for instance in:
DECLARE @Track geography;
SET @Track = dbo.MakeLineString(10)
SELECT @Track AS Shape, @Track.STNumPoints() as NumPoints;

Now I want to fill all fields Geo in the table Track with T-sql. In principle it's easy, my code is:
update Track
set Geo = dbo.MakeLineString(TrackId)

However, it does not run.
The error message is:
Msg 6522, Level 16, State 1, Line 1
A .NET Framework error occurred during execution of user-defined routine or aggregate "MakeLineString":
System.Runtime.InteropServices.COMException: The new transaction cannot be included in the list for the specified transaction coordinator.
System.Runtime.InteropServices.COMException:
   At System.Data.Common.NativeMethods.ITransactionJoin.JoinTransaction(Object punkTransactionCoord, Int32 isoLevel, Int32 isoFlags, IntPtr pOtherOptions)
   At System.Data.OleDb.OleDbConnectionInternal.EnlistTransactionInternal(Transaction transaction)
   At System.Data.OleDb.OleDbConnection.Open()
   At MyClass.MakeLineString(SqlInt16 TrackId)

I have no idea what's wrong.


texts: Split shapefile features and transpose attributes

texts:
QGIS v3.4

I have a single shapefile layer with overlapping features. I wish to intersect the features against themselves and transpose the summary of attributes in the output. The image below will help explain the goal.

What sequence of processing tools should I use to do this? The output should leave me without any overlapping features.


texts:
I’m struggling with returning data from the following query. It should be returning all data within a 5k radius of a given Longitude, Latitude point.

SELECT stuff
    FROM MyGeoTable
    WHERE ST_DWithin(ST_GeographyFromText(‘SRID=4326;POINT(LONG,LAT)’), geog,5000)

It is not returning all data it should be best guess there is an issue with the Geography column  this was created from a geometry column (type Geometry) and is a collection of lines and polygons all on WGS84.

I tried to recreate the Geography column thus

UPDATE MyGeoTable SET geog = (ST_Transform(geom,4326));

This gives

  ERROR: Geometry type (LINSTRING) dose not match column type (Polygon)

My question is twofold

how can recreate the geography column, and check it has worked on every row?
If this is not possible How might perform a similar query just using the geometry the 5k radius doesn't have to be exact but is from a user input!


texts: Tiled WMS server producing no data on Openlayers map?

texts:
I've been investigating setting up a server with OpenLayers which will display tiled data from a datacube-ows wms server. At the moment I have an index.html file which initializes the server and a wms.js javascript file which sets up the map I'd like to display.

It doesn't display any errors currently, but when I go to localhost:1234, there is an empty white page with no data on it.

I'm pretty new to wms, openlayers and web visualisation in general so I imagine there is probably something relatively simple I'm missing, but I haven't been able to find any answers from Googling so far. I reworked the wms.js file in the style of GeoServer WMS layer won't load on OpenLayers 5 but it had no effect.

index.html:

<!DOCTYPE html>

    <meta charset="utf-8">
    <title>Test Map</title>

wms.js:

import Map from 'ol/Map.js';
import View from 'ol/View.js';
import TileLayer from 'ol/layer/Tile.js';
import OSM from 'ol/source/OSM.js';
import TileWMS from 'ol/source/TileWMS.js';

var layers = [
  new TileLayer({
    source: new OSM()
  }),
  new TileLayer({
    source: new TileWMS({
      url: 'http://localhost/datacube_wms',
      params: {'LAYERS': '10DayIndices:ndvi', 'TILED': true},
      serverType: 'datacube_wms',
      projection: 'EPSG:4326'
    })
  })
];

var view = new View({
  projection: 'EPSG:4326',
  center: [48.031532, 102.376126],
  zoom: 2
});

var map = new Map({
  target: 'map',
  layer: layers,
  view: view
});

The WMS does work and below is an example of a successful GetMap request;

http://localhost:5000/?service=WMS
&request=GetMap
&version=1.3.0
&layers=10DayIndices
&styles=ndvi
&crs=EPSG:32648
&bbox=399600,5299260,499560,5399220
&width=4998
&height=4998
&format=image/png


texts:
I have a problem with turning on or off ArcMap editor, using arcpy.

When I try to turn on editor on arcpy I use this code on ArcMap Python window

workspace = r'C:\tests\test.mdb'
edit = arcpy.da.Editor(workspace)
edit.startEditing()
edit.startOperation()

But it doesn't activate editor. When I do this process with using Editor - Start Editing on Arcmap, I try to test it with this code :

workspace = r'C:\tests\test.mdb'
edit = arcpy.da.Editor(workspace)
xa = edit.isEditing()
print xa

It returns False, Although editing is on.

I need to deactivate editor in my Python add-in script. So, i have to use arcpy.


texts: Raster - Error NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

texts:
I am brand new to using geographic data in R and I keep running into an error:

  NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

I believe that this is due to an old version of PROJ (4.9.0), and I have found a good deal of discussion surrounding this topic but I have been unsuccessful in fixing this issue - mostly due to that I don't fully understand the answers posted on this topic i.e. "add missing file to proj folder" What file? Is my project folder the same as my working directory?.  I tried reinstalling rgdal, with the following commands:

install.packages("rgdal",
                 configure.args = c("--with-project-include=/p/bin/proj4/include",
                                    "--with-proj-lib=/p/bin/proj4/lib",
                                    "--with-proj-share=/p/bin/proj4/share/proj"))

But I get the same error.
I am a Mac user. Can someone walk me through how to deal with this error?


texts: I have problem where the street text elements overlapping place text elements in QGIS. Both elements are in separate layers. How to express in QGIS 3 that street name layer wont overlapping place name layer? Do you have idea or template about advanced street name labeling in QGIS?


texts: Python error when using r.reclass on raster with no CRS

texts:
I have a set of GRID data (directory of .adf files) that I'm trying to reclassify in QGIS.

I'm using the GRASS r.reclass function to try to reclassify the layers, but GRASS keeps failing and throwing a python error. When looking through the log, it seems to be an issue with the projection. The data I have "reports no coordinate system" and thus defaults to WGS-84 when imported to QGIS. Looking through the log, the first error it gives is this:

GRASS_INFO_WARNING(15749,1): Over-riding projection check
GRASS_INFO_END(15749,1)

GRASS_INFO_ERROR(15749,2): Illegal latitude for North
GRASS_INFO_END(15749,1)

GRASS_INFO_ERROR(15750,1): Invalid input
GRASS_INFO_END(15750,1)

I recognize that this is just a partial section of code and may not be that helpful, but any ideas on why GRASS/QGIS is throwing this error? I was able to successfully reproject other layers (with CRSs), so I don't think it's a compatibility issue between GRASS/QGIS. I also tried converting the data directory of .adf files to a .tif with GDAL thinking that might help, but no luck.

If the full log is helpful, I can update this post with that.


texts: GRASS GIS: how to run g.run_command("v.extract") in Python with sqlite, smooth-spline in GRASS GIS?

texts:
I would like to solve several problems by using a python script in
GRASS GIS:

I already uploaded collected GPS-waypoints and tracks ("v.in.org"), patched several GPS-vector data ("v.patch") and updated the attribute tables ("db.in.org","v.db.join"). I want generate a trail network, every trail is seperated in different sections (point vector data). In between the sections I tracked myself every minute (point data).

Now, I want to extract the waypoints of one trail (new vector) and connect the sections. Therefore, I tried:

>>>sqlite="""substr(cmt,1,4)='R001'"""
>>>print sqlite
>>>substr(cmt,1,4)='R001'
>>>g.run_command("v.extract", overwrite=True, input="rec_RG", output="test", where=sqlite)

0

As a second step I would like to connect the point data. I used the command "v.to.lines", but somehow GRASS integrates the area in between all points. How do I only connect the point-vector data from th beginning to the end? What do I need to add? I could not find any specific information in the manuel.

I tried:

g.run_command("v.to.lines", overwrite=True, input="spl_"+Nr_RG[i], output="spl_line"+Nr_RG[i])


texts: How to properly use Open Street Map's API

texts:
I'm using GeoTiler in my Python project to display a map view of some data. GeoTiler achieves this by querying some OpenStreetMap API (I assume the Overpass API). I have no intention of ever publicizing my code, so GeoTiler will likely make no more than a hundred or so requests per day. My code works some of the time, while failing other times for no appearent reason. I suspect that this is due to the fact that I'm using GeoTiler without an API key.
My code:
import geotiler

map = geotiler.Map(center=(-6.069, 53.390), zoom=16, size=(512, 512))
image = geotiler.render_map(map)
image.save('map.png')

Do I need an API key in order to use GeoTiler this way? If so, how do I obtain one? Are there other terms, conditions or pitfalls I need to be aware of?


texts:
I have georeferenced an old jpg map against Bing Road maps (OpenLayers plugin) and Raster>Georeferencer. I'm trying to use EPSG:4326 because my PostGIS database uses that. My PostGIS info shows coordinates in decimal degrees where they should be. But my georeference jpg show coordinates as -13000000 (not all zeros, but 13 million something), 4000000 (4 million something). Seeing EPSG:4326 as the CRS. What are these large coordinates and what do they mean?

Main problem is that the georeferenced jpg isn't where it's supposed to be. I have no idea where it is. I'm using zoom to layer to see it.

Side question. Is is possible to zoom to two layers at once? Would help with problems like this.

I'm using latest QGIS on a Mac.


texts: Getting all of combinations of features found on two layers using QGIS

texts:
This is a problem that seems like it should be simple, but it's a little difficult to explain. I have two polygon layers, each with a few hundred features (specifically the districts of England and the constituencies of England). What I want is to produce a third polygon layer, which would take in all of the combinations of districts and constituencies.
See the example below for a much simpler version, the first layer has features 1 and 2, the second layer has features 3 and 4. In this case you would want to get four features out -the intersection of 1 from the first layer and 3 from the second, 1 from the first and 4 from the second, 2 from the first and 3 from the second, 2 from the first and 4 from the second.
Is this possible, and is there a more elegant way of phrasing this that I'm not seeing?


texts:
I have exported osm data using osm2pgrouting in the database, it has SRID 4326. I want to change projection of all the tables created to a particular SRID in UTM.
I read the post FIrst and Second, but these are only changing SRID withour transforming the coordinates. I am using the command

select UpdateGeometrySRID('Schema Name', 'mytable', 'the_geom', newSRID) ;


texts: LAS Point Statistics as Raster in ArcGIS Pro fails to write raster

texts:
In the past I was able to use the ArcGIS tool LAS Point Statistics as Raster, but it is failing to work in ArcGIS Pro 2.2 or ArcGIS 10.7 on Windows 10.

After creating an LAS dataset, issuing the following command results in an error message:

arcpy.management.LasPointStatsAsRaster('lasdataset', 'outraster', "POINT_COUNT", "CELL_SIZE", 10)

The ArcGIS error message is 000210, "Cannot create output". I have tried to write to both a geodatabase and a folder.

When I use the following command, the function works:

arcpy.management.LasPointStatsAsRaster('lasdataset', 'outraster', "POINT_COUNT", "OBSERVATIONS", 10)

I have tried it in ArcGIS 10.7, Pro 2.2 both in the GUI and using arcpy. I really need a cell size of 10 m, and the observations cell size is huge.


texts: Qgis2leaf doesn't show layers with scale-dependent visibility

texts: When I add layers to the QGIS2LEAF dialog box it accepts it and opens up the preview in my web browser, but when I starting looking for the layers, in this case points for campgrounds, I dont see any of them. I have them on a 1:145,000 Scale Dependent Visibility minimum and 1:1 Maxium, but when I zoom in to where they should appear I dont see anything. Feel Free to ask for more information, I will add a couple photos the way I have it set up


texts: Minidump on exit in QGIS (2.8.2- Wien) 64-bit version when using GDAL in plugin

texts:
I am working on a plugin for QGIS that is working fine on 32-bit versions but I am having some problems while using it on the 64-bit standalone installation package.
I have tried on several machines and on all of them QGIS locks up and creates a minidump on exit when certain parts of the code in the plugin has been ran. The plugin still works, the only problem seems to be the crash dump on exit.

As far as I have been able to tell this line of code using the gdal-module seems to be causing the crash dumps:

    gdal_dsm = gdal.Open(filepath_dsm)
    dsm = gdal_dsm.ReadAsArray().astype(numpy.float) - this line of code is the offending one.
    ...
    gdal_dsm = None
    dsm = None

Am I doing something wrong or is it a problem with the 64-bit standalone installation of QGIS? Is there anyway to avoid the crash dump on exit for potential users of the plugin that might be using this installation package?


texts: Leaflet Vector Layers (lvector.js) not working in layercontrol

texts:
Im using lvector.js from Jason Sanford to draw Postgis-data as geojson on top of a leaflet map.
With the following code, the layers 'layer_n2k_do' is drawn correctly on the map:

    var layer_n2k_do = new lvector.PRWSF({
    url: "http://map.domain.com/rest",
    geotable: "n2k_do",
    fields: "gid, sitecode, sitename, surfha",
    where: "",
    uniqueField: "gid",
    srid: 2169,
    showAll: false,
    scaleRange: [14, 20],
    map: map,
    singleInfoWindow:true,
    popupTemplate:"

{sitecode}
 {sitename}

{surfha} ha

",
    symbology: {
        type: "single",
        vectorOptions: {
            fillColor: "orchid",
            fillOpacity: 0.1,
            weight: 2.0,
            color: "orchid",
            opacity: 0.8,
            clickable: true
        }
    }
});

When I don't put the layer directly to the map, i can later draw it on the map with:

layer_n2k_dh.setMap(map);

I would like to use a layercontrol to collect all my overlays, so that I could toggle their visibility just by activating/deactivating them in the control.
I've tried the Leaflet.groupedlayercontrol, but the layers doesn't appear on the map?

var basemaps = {
        OSM: L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
    attribution: '&copy; OpenStreetMap contributors'
    }),
        'Google_Hybrid': L.tileLayer('http://{s}.google.com/vt/lyrs=s,h&x={x}&y={y}&z={z}',{
    maxZoom: 20, subdomains:['mt0','mt1','mt2','mt3']
    })
    };

    var groupedOverlays = {
      "Natura 2000": {
        "Habitat-Dir.": layer_n2k_dh,
        "Birds-Dir.": layer_n2k_do
      },
      "Borders": {
        "Municipalities": layer_mun
      }
    };

    L.control.groupedLayers(basemaps, groupedOverlays).addTo(map);
    basemaps.OSM.addTo(map);

The layers are listed in the control, but they aren't added to the map.

Does anyone have an idea how to get thar work? Or is there an other control that someone knows top work with lvector.js-layers?


texts:
What is the best aproach to make a QGIS plugin using PyQt (PyQGIS), and use that same app and integrate with a django project, in order to put that app on the web with the minimum replicated code possible.

I understand that this question may be a little bit vague, but what I intend to know is the possible approaches to developing web based applications using QT with python instead C++.

I have found this one possible aproach which is QtWebKit.


texts: In R's raster::terrain function, which way is north?

texts:
The terrain function in R's raster package calculates slope and aspect from an elevation raster:  see http://www.inside-r.org/packages/cran/raster/docs/terrain.

With the unit="degrees" option, the function returns aspect values in [0, 360].  Which way is north?  Which way is west?  I am not an arcgis user, but arcgis documentation says

  Aspect identifies the downslope direction of the maximum rate of
  change in value from each cell to its neighbors. It can be thought of
  as the slope direction. The values of each cell in the output raster
  indicate the compass direction that the surface faces at that
  location. It is measured clockwise in degrees from 0 (due north) to
  360 (again due north), coming full circle. Flat areas having no
  downslope direction are given a value of -1.

Is that standard?  Does raster::terrain in R follow the same convention, i.e. do 0 and 360 both mean north?

Snippet of R code:

> library(raster)
> swiss_elevation <- getData('alt', country='CHE')
> raster_brick <- terrain(swiss_elevation, opt=c('slope', 'aspect'), unit='degrees')
> raster_brick$elevation <- swiss_elevation
> plot(raster_brick)
> summary(raster_brick$aspect[])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
   0.00   88.58  179.60  183.60  283.70  360.00   81354

This little example suggests that 0 is indeed north, 90 is east, et cetera:

library(raster)

test <- raster(ncols=100, nrows=100)
test[] <- 1
max_i <- 10
for(i in seq(max_i, 1)) {
    test[seq(50-i, 50+i), seq(50-i, 50+i)] <- max_i - i + 1  # "Mountain" at center of raster
}
plot(test)
plot(terrain(test, opt="aspect", unit="degrees"))  # What causes the bands?
plot(terrain(test, opt="slope"))

Does any of the R raster documentation point that out, or is it so obvious that it doesn't need to be stated?


texts: How to improve GeoServer WCS performance?

texts:
I have a Geoserver installation which I am using to serve a large area of 2m LiDAR DSM as a Web Coverage Service (WCS). I'm finding the performance is pretty slow (at least when using the WCS in QGIS) and am looking for advice on tuning it or in case I've misconfigured something.

I'm not particularly familiar with WCS. The bbox of the data is about 100km by 70km. I expect the response to be very slow at small scales and thus have only been testing it when at high zoom levels, but it is still ~1 or ~2 mins a query and can't handle concurrent queries.

I've configured the data as an ImageMosaic store of ~3000 geotiffs granules. Each tif covers 1km and was configured with the -co "TILED=YES" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" parameters in GDAL when they were created. I just use the default generated index and properties file when I load these as an ImageMosaic in Geoserver.

Each granule is about 2 mb so I would've thought if calling getCoverage for a few tens of these shouldn't really lock up the server up very along.

The server is Ubuntu, Geoserver 2.5.1, 3gb ram,

The tomcat settings are:

JAVA_OPTS="-Djava.awt.headless=true -Xms1024m -Xmx1024m XX:PermSize=128m -XX:MaxPermSize=512m -XX:+UseConcMarkSweepGC"

Am I expecting too much? Is there anything else I could tweak to improve performance?


texts:
How can I read a gdal vrt file with GDALOpen function?

I have a vrt file built with gdalbuildvrt and I want to read it using GDALOpen  but It returns a null dataset, Is it even possible?
In this code poDataset is null

GDALDataset  *poDataset;
        GDALAllRegister();

           poDataset = (GDALDataset *) GDALOpen( "home/m/3857/srtm_merged.vrt", GA_ReadOnly );
        if( poDataset != NULL )
        {

}

It seems vrt files are ony supported in a few gdals tools.

Here is VRT:

<VRTDataset rasterXSize="21746" rasterYSize="17369">
  <SRS>PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</SRS>
  <GeoTransform>  5.0093307430080175e+06,  5.1196462744776348e+01,  0.0000000000000000e+00,  4.3927264476078171e+06,  0.0000000000000000e+00, -5.1196462744776348e+01</GeoTransform>
  <VRTRasterBand dataType="Int32" band="1">
    <NoDataValue>-9999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="5344" RasterYSize="6593" DataType="Int32" BlockXSize="5344" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="5344" ySize="6593" />
      <DstRect xOff="0" yOff="0" xSize="3624.28574350654" ySize="4471.35402450199" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">2.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="5529" RasterYSize="6439" DataType="Int32" BlockXSize="5529" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="5529" ySize="6439" />
      <DstRect xOff="0" yOff="13147.5894347461" xSize="3624.84957383552" ySize="4221.45169215535" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">3.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="5470" RasterYSize="6489" DataType="Int32" BlockXSize="5470" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="5470" ySize="6489" />
      <DstRect xOff="10871.7953859191" yOff="4468.83771045391" xSize="10874.354045833" ySize="12900.1249366381" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">4.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="5344" RasterYSize="6593" DataType="Int32" BlockXSize="5344" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="5344" ySize="6593" />
      <DstRect xOff="10871.7953859191" yOff="0" xSize="3624.28574350655" ySize="4471.354024502" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>

well to others may face the same problem, there must be a / at the beginning of vrt path

GDALDataset  *poDataset;
        GDALAllRegister();

           poDataset = (GDALDataset *) GDALOpen( "/home/m/3857/srtm_merged.vrt", GA_ReadOnly );
        if( poDataset != NULL )
        {

}


texts: Building version control system for website, like the one used on Wikipedia but for geospatial data?

texts:
I'm working on a website permitting to display and modify geospatial data. So far, I'm using this solution:

My data is stored in PostgreSQL with the extension PostGIS
It is posted in GeoServer
My website uses OpenLayers to request the tiles I need from GeoServer with WMS and WFS

Now I would like to keep the memory of the previous versions of the geospatial data. On Wikipedia, users can read an article, modify it or check the discussion page of the article to see the previous commits. I would like to have the same possibilities with the geospatial data:

A page to see an element (let's say a multipolygon)
A page to modify it and a page to see all the commits made with this element (with the author, the date and the size of each commit
And the possibility to click on a commit to see the geospatial data of the element at the time of the commit)

I tried to see what existed on the internet and it sounds like what I need is a version control system (like the one on Github) for geospatial data. Geogig (http://geogig.org/) seems to be the one. But I'm not quite sure I understand all this well. According to my needs, do you think a version control system and more especially Geogig could be a solution?
Also, even if Geogig can do the job, I'm not sure How I could make it work. As I understand it, I would put an extension to use Geogig in GeoServer and Geogig would be an intermediate between PostgreSQL and GeoServer. Does it mean I would keep the connection PostgreSQL<->GeoServer for the current data (so for the display page and the modification page of an element) and the connection PostgreSQL<->Geogig<->GeoServer for the past commits (so for the discussion page of an element)?


texts:
I've written a script in a python toolbox. That script creates a folder in the user defined workspace for all the output files that are created during the script process. When the script is done, I don't need these files anymore because I've saved my results in an extra folder. I want to delete this temporary folder and all of its files and wanted to use the delete_management tool like this:

# Set Workspace
    arcpy.env.workspace = parameters[2].valueAsText

    # Set Temporary Output
    arcpy.CreateFolder_management(arcpy.env.workspace, "Temp")

    # Delete temporary Files
    arcpy.Delete_management("Temp")

But that clearly doesn't work and it is not even giving me an error message. The delete line just has no effect. What am i doing wrong here?

edit: I'm using ArcGIS 10.1 and the parameters[2] string is the users input.


texts: For a non-geocentric ellipsoid, how does the Lat/Long/Ell Ht & semi-major axis & flattening implicitly define the centre of the ellipsoid?

texts:
For example if we compare the AGD84 (Australian Geodetic Datum) with GDA94 (Geocentric Datum Australia), the centre of the AGD Datum is 200 metres different to GDA, because it is a best fit to the geoid.

But how is this difference realised (implicitly) given only the Lat/Long/ Ellipsoidal Height of a Point (Johnstone), the semi-major axis (a) and the flattening (f)? I realise that (a) is different by 23 metres, but it is the angle between the 2 datums that makes the most difference, this is what I am trying to get a handle on.

Asking another way, how could I calc the ECEF coords for the centre of the AGD, to see how different it is to GDA?

Asking yet another way, how to determine the shift & angles between the axes of of AGD & GDA?

I already have the GDA manual - http://www.icsm.gov.au/gda/gdatm/gdav2.3.pdf
 and am working through writing C++ code for this. I also understand how AGD is not very consistent, so a conversion between the 2 datums requires the use of a model, it is just that I am trying to truly understand the mathematical difference between the two.

Sorry for the particularly nerdy question, & thanks in advance for any answers.


texts: Creating coloured map from Excel using ArcGIS for Desktop?

texts:
I have an Excel file with some data and I want to create a map and present data with colours. I do everything well except for I present my data with dots. Instead, I need to cover the free space with an average of neighbour dots so as to have a solid map. Check out the photos to understand what I am saying.

This is what I have:

and this is what I want to make:

All the above has to be done in ArcMap.


texts: Jar conflicts between tomcat webapp makes geoserver doesn't work

texts:
I've been developing a simple app that uses geoserver inside tomcat, I was intending to register data to PostGis using java code and retrieve on Openlayers using geoserver. My webapp use a Postgres (https://jdbc.postgresql.org/download.html) and the PostGis (https://github.com/postgis/postgis-java/tree/master/jdbc) driver (both in tomcat lib directory), but Geoserver wasn't working properly, after looking in internet I found this:https://stackoverflow.com/questions/6981564/why-must-the-jdbc-driver-be-put-in-tomcat-home-lib-folder.

I suspect there is some jar conflict between the tomcat and geoserver probably with the PostGis driver and maybe geotools jar. With this conflict I'm not able to create DataStores in geoserver with my PostGis database. If I remove the PostGis jar from tomcat lib directory, GeoServer works correctly.
I'd like to know if the problem is with Postgis driver, and if is possible to use them together somehow in the same tomcat instance.

[UPDATE]
I'm using Tomcat 7.0, geoserver WAR version 2.10. Postgresql jdbc driver v9.4.1211 (the same in geoserver to avoid conflicts) and
PostGis jdbc driver 2.2.2, I followed as explained in github repository Readme.md: "You can use the org.postgis. DriverWrapper as replacement for the
  jdbc driver. This class wraps the PostGreSQL Driver to
  transparently add the PostGIS Object Classes.".

When I try to create a new DataStore inside geoserver to my PostGis database a red message shows up just over the button "Save"

Error creating data store, check the parameters. Error message: Unable to obtain connection: Cannot create PoolableConnectionFactory (
Something unusual has occurred to cause the driver to fail ...)

Inside eclipse, the error log is like this:

  [data.store] - Error obtaining new data store java.io.IOException     at
  org.geoserver.catalog.ResourcePool.getDataStore(ResourcePool.java:676)
    at
  org.geoserver.catalog.impl.DataStoreInfoImpl.getDataStore(DataStoreInfoImpl.java:38)
    at
  org.geoserver.web.data.store.DataAccessNewPage.onSaveDataStore(DataAccessNewPage.java:95)
    at
  org.geoserver.web.data.store.AbstractDataAccessPage$1.onSubmit(AbstractDataAccessPage.java:165)
    at
  org.apache.wicket.ajax.markup.html.form.AjaxSubmitLink$1.onSubmit(AjaxSubmitLink.java:110)
    at
  org.apache.wicket.ajax.form.AjaxFormSubmitBehavior$AjaxFormSubmitter.onSubmit(AjaxFormSubmitBehavior.java:215)
    at
  org.apache.wicket.markup.html.form.Form.delegateSubmit(Form.java:1307)
    at org.apache.wicket.markup.html.form.Form.process(Form.java:974)

I'd like to post printscreens for the errors but I'm not able because of my reputation.
[ENDUPDATE]


texts: Can a map be added to a report generated through ExportReport using Python?

texts: I am currently able to generate a PDF report using arcpy's ExportReport that contains attribute information for records within an attribute table, but is there a way to add a map to that report, as another page, using python? I can't seem to find an documentation on this.


texts: Is it possible to use assign_vertex_id() POLYGONs?

texts:
I'm trying to use pgrouting in order to create a graphe and I have two questions :

Is it possible to use the function assign_vertex_id() with a POLYGON geometry type? I always have the following error when I try to use this fonction :

psql (9.0.3)
Type "help" for help.

bpi=# SELECT assign_vertex_id('rangee', 0.00001, 'geometrie', 'id');
NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"

CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"

PL/pgSQL function "assign_vertex_id" line 14 at EXECUTE statement

ERROR:  query string argument of EXECUTE is null
CONTEXT:  PL/pgSQL function "assign_vertex_id" line 31 at EXECUTE statement

the table rangee is :

bpi=# \d rangee

          Table "public.rangee"

    Column    |         Type          | Modifiers
 -------------+-----------------------+-----------
  id          | integer               | not null
  zone_id     | integer               |
  geometrie   | geometry              |
  description | character varying(50) |
  source      | integer               |
  target      | integer               |

Indexes:
 "rangee_pkey" PRIMARY KEY, btree (id)
 "rangee_spatial" gist (geometrie)

Triggers:

rangee_geom_srid BEFORE INSERT OR UPDATE ON rangee FOR EACH ROW EXECUTE PROCEDURE geometrie_srid()

I don't understand this mistake. Could you help me please?

My second question is about pgrouting. Is it possible to create a dual graph with this tool?

Thanks a lot  and I'm sorry for my english.


texts: How to import multiple TIFF files via r.in.gdal in GRASS GIS?

texts:
I have a 18-year time series of the MODIS with temperature data, however there are more than 750 TIFF files.

How can I do to import all at once? Importing the files one by one is a very difficult task.

How do I do this using the function r.in.gdal in GRASS GIS?


texts: Is there an alternative function which does the same thing as clip raster by masked layer on QGIS?

texts: I am having some technical problems with the clip raster by masked layer and I faintly remember using a different method a long time ago for a similar problem. Anybody familiar with another function?


texts:
I'm trying to label visited countries in a world map shapefile, what is wrong with my rule :

if( "Visited" = 'yes', "NAME", NULL)


texts: Visibility On-Off of Label Layer in arcgis javascript 3.7

texts:
I am using ArcGIS JavaScript 3.7 API's and i am using agsjs.dijit.TOC control and it is working fine but when i added Label Layer in my map that time Label layer is coming fine but TOC is not working Properly.

While Pressing the checkbox to On-Off the Layer it is working fine but regarding the feature layer i added one Label layer that is not getting On-Off

Can anyone tell me how to Switch-Off or On the Label Layer.

Thanks in Advance.

I am Pasting some code..

    Var load = function () {
    require(["dojo/_base/Color", "esri/symbols/TextSymbol", "esri/renderers/SimpleRenderer", "esri/renderers/ScaleDependentRenderer"],
    function (Color, TextSymbol, SimpleRenderer, ScaleDependentRenderer) {
    var statesColor = new Color("#666");
    var fieldLabel = new TextSymbol().setColor(statesColor);
    fieldLabel.font.setSize("4pt");
    fieldLabel.font.setFamily("arial");

    statesLabelRenderer = new SimpleRenderer(fieldLabel);
    return statesLabelRenderer;
  });
    }
    var layers = [];
    var layer = new esri.layers.FeatureLayer(url, {
        mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
        outFields: ["*"]
    });
    layers.push(layer);

    layerInfo.push({ layer: layer, title: item["LookupDisplayDesc"], boxLabel:     item["LookupDisplayDesc"], id: layer.id, checked: layer.visibleAtMapScale, slider: true });

if (item["LookupLongDesc"] != null && item["LookupLongDesc"] != "") {
    var fieldNames = item["LookupLongDesc"].split(",");
    var labelFields = "";
for (var i = 0 ; i < fieldNames.length ; i++) {
    labelFields = labelFields + "${" + fieldNames[i].trim() + "},";
}
    labelFields = labelFields.slice(0, -1).trim();
    load();
    var labelLayer = new esri.layers.LabelLayer();
    labelLayer.addFeatureLayer(layer, statesLabelRenderer, labelFields);
    layers.push(labelLayer);
    }
}
map.addLayers(layers);

 dojo.connect(map, 'onLayersAddResult', function (results) {
                if (layerInfo.length > 0) {
                    var legendDijit = new esri.dijit.Legend({
                        map: map,
                        layerInfos: layerInfo
                    }, "legendDiv");
                    legendDijit.startup();
                }
                toc = new agsjs.dijit.TOC({
                    map: map,
                    layerInfos: layerInfo
                }, 'LayerDiv');
                toc.startup();
            });


texts: Identifying polygons overlapping several administrative levels

texts:
I have information about grid population in a country assigned to polygons; the borders of these polygons are defined by night light intensity. In this way I am able to identify the "real" extent of different urban agglomerations that incidentally go beyond the official administrative borders.

As an example, consider the screenshot attached:

The yellow polygon represent an urban cluster, the black lines are the administrative border of ADMIN 3; the red lines are the borders of ADMIN 2 (where ADMIN 2>ADMIN 3 in the governance hierarchy). Borders from ADMIN 3 and ADMIN 2 come from two different shape files.

Now, I would like to identify with which ADMIN 2 and/or ADMIN 3 this polygon (and similar others in the country) overlaps. In other words, I would like to create an attribute table that tells me:

polygon with FID 1 is common to/overlaps with ADMIN 2 "xy" and ADMIN 3 "zw".

Is it possible to obtain something similar in QGIS 2.14.3 for Mac?


texts: Calculate the majority of one polygon in another (fishnet)

texts:
I am looking for a method to calculate the majority of one polygon file (a map of ecoregions) in a sampling grid (fishnet) in ArcGIS.  The purpose is to assign each sampling grid (fishnet) to an ecoregion.  The determining factor will be the majority of area.

I had a decent method of working out the problem using Tabulate Area tool, but the numbers were almost all identical indicating something weird was going on. I didn't trust it.  I was also playing around with the intersect tool then dissolve, but couldn't find a way keep the attribute data I needed (in this case, the ecoregion I want to label).

Problem: Label each sampling grid by the majority area of ecoregion within the sampling grid.

Can someone help me out?


texts: Filter two image collections by the same aqusition date in Google Earth Engine

texts:
I am trying to compare images aquired from two different sensors: MSI on board of Sentinel-2 and MODIS on board of Terra, using Google Earth Engine.
MODIS provides his own band of the parameter I am instered in at 500 km resolution  concentrantion meanwhile I created my own algorithm in Sentinel-2 at finer resolution.
I want to compare the results over a specific area but first I want to filter the two collections by selecting only images that were aquired approximately on the same date, it only matters the YYYY-MM-DD.
Specifically TERRA has a daily temporal resolution over the selected area while Sentinel-2 5 days. Basically I want to reduce the MODIS collection selecting only those images overlapping temporally with Sentinel-2.
Note that I will nead to map the results.
I don't have any idea on how to do it, I came across the function 'inner.Join' but I don't know how to use it.
I am not planning to merge them into one Image collection, I want to keep them separetely


texts:
I am having trouble displaying popups on points plotted on my OpenLayers 3 map from two different WFS sources. Both are PostGIS tables coming from GeoServer.

var popup = new ol.Overlay.Popup();
map.addOverlay(popup);
map.on('click', function(evt) {
    // Hide existing popup and reset it's offset
    popup.hide();
    popup.setOffset([0, 0]);

    // Attempt to find a marker from the planningAppsLayer
    var feature = map.forEachFeatureAtPixel(evt.pixel, function(feature, layer) {
        return feature;
    });
    if (feature) {
        var coord = feature.getGeometry().getCoordinates();
        var props = feature.getProperties();
        var info = "
" + props.addr + "
";
            info += "
" + coord + "
";
            info += "
Status: " + props.addr + "   " + props.stuff + "
";
        // Offset the popup so it points at the middle of the marker not the tip
        popup.setOffset([0, -22]);
        popup.show(coord, info);

    }

Addr comes from the layer fpal, and stuff comes from the layer planningAppsLayer. I was wondering if there was a way to create popups with completely different popup for the information from each layer. Not sure if I am making sense.


texts: Reading output spatial reference from geoprocessing service in Python toolbox (.pyt) service?

texts:
I know I could just create a custom parameter but it seems there's already a built in method here to specify the output spatial reference - I just have no idea how to read it.

I've tried checking arcpy.env.outputCoordinateSystem but it's null at runtime.


texts: With ArcGIS Server's REST API, can I request a non-default symbol size when requesting a legend JSON?

texts:
I am using the ArcGIS Server REST API to retrieve legend information for one of our map services. I am requesting that information as a JSON object, so my HTTP request URL looks like this:

http://<host>/ArcGIS/rest/services/<mapServiceName>/MapServer/legend?f=json

The legend symbol images I get back from this are all 20×20 pixels in size (which seems to be the default). Guessing from the relevant bit of Esri's documentation, there is no way to request a different symbol image size (e.g. larger images for high-resolution printing). That would, however, be possible using the SOAP API (I believe by specifying MapServerLegendPatch objects).

Is there any way to request a different legend symbol image size using the REST API?


texts: Why image do not color match when many Landsat images are mosaicked in ERDAS?

texts:
When I try to mosaic Landsat images in ERDAS I get as result image with discordance in colors.

I tried histogram matching but no success.

The images were captured between March and May.


texts:
I'm trying to delete a map that I used as a mask before, but I cannot delete it using g.remove. It displays this:

g.remove rast=raster_Mnorzagaray@PERMANENT
[raster_Mnorzagaray@PERMANENT@PERMANENT] is a base map. Remove reclassed map first: MASK@PERMANENT


texts: Adding layer to group in layers panel using PyQGIS?

texts:
I am trying to add a group to the layer panel, then place a vector layer (vectorLayer) into this group from a python script. The code I am using is:

groupName="test group"
root = QgsProject.instance().layerTreeRoot()
group = root.addGroup(groupName)

group.insertChildNode(-1, QgsLayerTreeLayer(vectorLayer))
QgsMapLayerRegistry.instance().addMapLayer(vectorLayer,True)

My issue is that this code always adds two copies of the vector layer. One of which will be in the group, the other is outside the group.

However, when I change the code to the following with the intent of not adding a second layer, the layer inside the group becomes '(?)'.

groupName="test group"
root = QgsProject.instance().layerTreeRoot()
group = root.addGroup(groupName)

group.insertChildNode(-1, QgsLayerTreeLayer(vectorLayer))
# QgsMapLayerRegistry.instance().addMapLayer(vectorLayer,True)

Does this mean the layer inside the group is just a pointer to the layer outside the group?

How do I add a layer to a group without the duplicate layer?


texts: Quering out multiple values in MapInfo SQL Query

texts:
I am trying to filter out around 1000 rows and delete them in MapInfo as they are no longer required in the data set. If there was a smaller quantity then I would use an SQL Query but it has its limits and I can't get it to work.

The numbers/rows I need to filter out are currently in a Notepad in the following format 46061871, 46062951, 46063180, etc. and none of the other fields in the rows would allow me to delete them so it has to be done via the unique values.

Can anyone advise a workaround or solution?


texts: Creating dimension feature class that displays inches in ArcGIS Pro

texts:
I have a line feature class in State Plane survey feet. When I create a dimension feature class I simply cannot figure out how to get the dimension to display inches. There seems to be a button for inputting an expression (DimLenght*12 maybe), but it is greyed out. There is an option to set the display output to inches, but that doesn't seem to affect the final output - the dimension lines still show feet, but with an " suffix (I set this option).
Where is the option to change what is displayed on the dimension line?


texts:
I am trying to export our Property layer to DXF but am running into issues with labels.
I have tried each of the three ways that I have seen described about DXF export, and none of them work for me.

My issues are

a) Not all the labels are exported with the DXF and

b) Then when I open the file in AutoCAD, the font size is huge.

Neither my AutoCAD expert colleague nor I have been able to make the font size smaller in AutoCAD.

I know the documentation with the DXF export says that labels are not included in the export, but that was some time ago and I have hopes that this situation has changed.

Any suggestions?


texts: How to exclude several values in a definition query

texts:
I need to set up the definition query for a layer to exclude several types of values for a particular field. It is a status field that that has values such as Active, Inactive, Void, Received, etc. The definition query needs to be something along the lines of Status IS NOT "Inactive" OR "Received" OR "Void".

I've tried setting it up to Status IS NOT "Inactive" AND Status IS NOT "Received" AND Status IS NOT "Void"
and I've tried formatting it to IS NOT... OR...OR


texts: Automating calculation of NDVI with rasterio?

texts:
I want to automate this code below for other 99 images, the problem that I have is that each image is in its specific folder as it's mentioned in the picture below.

the images are very big and I want to clip them in a specific area.

Is that possible to make it automatic also ?

    import rasterio
import os
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np

os.chdir(r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP')

band4 = rasterio.open(r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP\clip_RT_L1C_T29SPR_A003995_20160328T111353_B04.tif') #red
band5 = rasterio.open(r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP\clip_RT_L1C_T29SPR_A003995_20160328T111353_B08.tif') #nir

plot.show(band5)
band4.dtypes

red = band4.read(1).astype('float32')
nir = band5.read(1).astype('float32')

ndvi=np.where(
    (nir+red)==0.,
    0,
    (nir-red)/(nir+red))
ndvi[:5,:5]

plot.show(ndvi)
plt.colorbar (ndvi) #I have also a problem with colorbar as it didn't show properly


texts: Creating join based on multiple fields using ArcGIS Desktop?

texts:
I have two feature classes in a File Geodatabase that I would like to join based on multiple fields.  I've searched this site and Google and all I have come up with was to use the Make Query Table tool.  I've tried this, but I keep getting an SQL error.  My SQL is pretty poor and I'm pretty sure I'm missing something.

I am aware I can create a new field and concatenate the values from my fields, but I would like to avoid this, if possible.

I'm using something that looks like this:

(Table1.Field1 = Table2.Field1) AND (Table1.Field2 = Table2.Field2) AND (Table1.Field3 = Table2.Field3)

When I verify the query, I get an error that says:

There was an error with the expression.
An Invalid SQL statement was used.
An invalid SQL statement was used. [Table1]

Also, if someone has another solution that doesn't use this tool, I'm happy to hear about it.


texts: Exporting to PDF using ArcPy of ArcGIS Pro?

texts:
I have a weird issue in exporting the layout to PDF
My code is really simple and I have tried in ArcGIS Python IDLE but I have came across the problem.

My Code and error is :

Code:

aprx = arcpy.mp.ArcGISProject("CURRENT")
Layout = aprx.listLayouts()[0]
Layout.exportToPDF(r"C:\1.pdf")

Error:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\utils.py", line 191, in fn_
    return fn(*args, **kw)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 975, in exportToPDF
    return convertArcObjectToPythonObject(self._arc_object.exportToPDF(*gp_fixargs((out_pdf, resolution, image_quality, compress_vector_graphics, image_compression, embed_fonts, layers_attributes, georef_info, jpeg_compression_quality, clip_to_elements), True)))
OSError: C\1.pdf


texts:
Our server is using UTM internally. We then send a large number of UTM Coordinates to display them in a custom made leaflet layer.

First we transform (unproject) them into "normal" geographical coordinates using Proj4Leaflet . Then we use leaflets latLngToLayerPoint() function in order to project them for our custom layer.

The latLngToLayerPoint() function is projecting the geographic coordinates back into point coordinates.

Is there a way to transform our UTM coordinates directly into points without making the detour by geographical coordinates? I see that there is a transform() function which is used by Leaflet internally to project "normal" latlng coordinates.

How can I get the necessary parameters to perform a transformation from EPSG 32633 to the internal 3857 projection?


texts: Watershed layers with names for each one for Europe or Sweden?

texts:
I am working on freshwater species and I would like, make a list of each watersheds used for each species.

So I am looking for a watershed maps which contain name of the watershed.

I found some good layers, but, they don't get names for each basins, only for the bigger ones.

I would like to know if someone can give me advices, links or good sources?

I need mainly the cover for Sweden at least, or Europe.

I searched this site, but I found most of the elements to create watersheds.


texts: Extracting road and building features from raster data using ArcGIS Desktop?

texts:
I work with raster data and my task is to digitize a big image. This consumes much time just for digitizing roads and buildings!

I'm looking for tools to simplify working on raster data to digitize features, such as automate road extraction, smooth features, etc.
For examples, check these videos: RoadTracker & Overwatch.

I have ArcGIS 9.3 and 10 but other suggestions are welcome too.


texts:
In reference to my previous question.

I am looking for the benefits of telling SQL Server your data is spatial.  For instance, I have a duplicate of that table generated by the Shape2SQL tool (see previous question).  The major difference is it is not a spatial table, but I have a field where I am storing the well-known text in lat/long (e.g. POINT(-82.2990401 38.9186892)).  I found I was able to run a similar query that worked to get me the data I was requesting.

declare @poly geometry = 'POLYGON((-82.8113 39.690, -81.145 39.725, -80.717 38.74, -82.813 38.614, -82.8113 39.690))'

SELECT * FROM [TABLE]
WHERE @poly.STIntersects(WKT) <> 0;

I am looking to see the real benefits to telling SQL Server that you're dealing with spatial data.  The only two things I can think of off the top of my head are spatial indexing and being able to visualize your data in the query results window within SQL Server Management Studio.

I am looking for the pros of SQL Server serving spatial data.  This is for a personal website, and I have no desire to invest in running a spatial server instance like ArcGIS Server or GeoServer.  I am curious if I will be better off leaving things in a table and just using WKT in my queries or if I should really be using a geom/geog object.


texts: How to update attribute table layout in QGIS Print Composer?

texts:
I have selected some columns from attribute table in print composer linked to layer A.

But when I open my layout in another project with that same layer A, it doesn't refresh/update

and I have to select columns again. Is there a solution to avoid this operation?


texts:
I am trying to add WMS layer on QGIS. It has been added but it says 'Protocol "" is unknown'.

The URL I am using are:
https://data.opendevelopmentmekong.net/geoserver/ODCambodia/Forest_cover_2009/wms?service=WMS&request=GetCapabilities&layers=ODCambodia%3AForest_cover_2009

https://data.opendevelopmentmekong.net/geoserver/ODCambodia/wms?service=WMS&version=1.3.0&request=GetCapabilities

How can I display it?


texts: Create Geospatial PDF (GeoPDF) is greyed out in PDF export options in QGIS 3.10

texts: I installed QGIS 3.10 standalone version today and I know it has a new feature introduced in this version to export GeoPDF. But for some reason, it is greyed out.


texts: Geomesa accumulo storing null Point for geometry

texts:
I am using Point geometry type in geomesa for strong lat-long details in feature type. In one of the use case while storing a new feature I don't have value for lat-long so I tried sending null object for createPoint function, but the library 'jts' throwing exception while processing null Geometry type. Here is the source code. By seeing the implementation of jtslibrary, It seems difficult for storing null in geomesa Point. What are the possibilities I can go for?
Here is the snapshot of my code:

sf.setAttribute(Attributes.GEOMETRY.getValue(), JTSFactoryFinder.getGeometryFactory().createPoint((Coordinate)null/*new Coordinate(event.lat, event.lng)*/));

Attaching logs for exception:

java.lang.IllegalStateException: getX called on empty Point
    at com.vividsolutions.jts.geom.Point.getX(Point.java:124)
    at org.locationtech.geomesa.utils.uuid.Z3UuidGenerator$.createUuid(Z3FeatureIdGenerator.scala:93)
    at org.locationtech.geomesa.utils.uuid.Z3UuidGenerator$.createUuid(Z3FeatureIdGenerator.scala:70)
    at org.locationtech.geomesa.utils.uuid.Z3FeatureIdGenerator.createId(Z3FeatureIdGenerator.scala:32)
    at org.locationtech.geomesa.index.geotools.GeoMesaFeatureWriter$.featureWithFid(GeoMesaFeatureWriter.scala:64)
    at org.locationtech.geomesa.index.geotools.GeoMesaFeatureWriter.writeFeature(GeoMesaFeatureWriter.scala:128)
    at org.locationtech.geomesa.index.geotools.GeoMesaAppendFeatureWriter$class.write(GeoMesaFeatureWriter.scala:194)
    at org.locationtech.geomesa.accumulo.data.AccumuloAppendFeatureWriter.write(AccumuloFeatureWriter.scala:19)
    at com.hps.GeomesaClient.insert(GeomesaClient.java:129)


texts:
I'm currently using the layer_property(layerName,'source') to produce a list of copyright statements for the Print Composer by putting that in an if() statement.

However, that isn't very flexible as it would need hardcoding for every layer and it wouldn't have checks to see if the copyright had already been printed for a previous layer.

The obvious solution is to put this into the function editor, but I can't find anything that would work in python/ function editor to pull the file path out.


texts: Creating and calling virtual environment to run ArcPy script?

texts:
I'm getting a new computer and am trying to set up my virtual environments the right way. In the past, I've just been using PyCharm and running everything through "C:\Python27\ArcGIS10.6" so it's a miracle I haven't broken ArcPy yet on my old computer.

I'm setting up virtual environments for:

just ArcPy with ArcMap (so, Python 2.7); and
Python 3.x for GDAL and other things

I'm following @CurtisPrice’s answer to Making separate Python installation that can call ArcPy? to set up virtual environments.

There is a simple example below, where I am writing and running a .py script in PyCharm.

Do I now just call @CurtisPrice's *.py file to set the virtual environment?

usercustomize.py   ## is this how I set the virtual environment???
import arcpy

myshp = r'C:\Users\name\path\myfile.shp'
arcpy.AddField_management(myshp, "new_field", "SHORT")


texts: Population calculation per grid using QGIS?

texts:
I have a ward map of a city along with the population in each ward.  I have divided the map into 1 km* 1 km grids.

How do I calculate the population per grid, taking into account that there can be more than 1 ward in each grid?


texts: PyQGIS Print Composer not exporting raster layer successfully (CRS in correct format)

texts:
I have a Python script that runs and creates statistics using zonal statistics sum tool. The python script adds each layer to the map successfully:
# Add Base Map
OSMBase = "....\\OSM_Base.xml"
rlayer2 = QgsRasterLayer(OSMBase, "OSM Base Map")
if not rlayer2.isValid():
    print("Layer failed to load!")
QgsProject.instance().addMapLayer(rlayer2)

# Add Population Raster
path_to_tif = "....\\LAC_Population.tif"
rlayer = QgsRasterLayer(path_to_tif, "LAC_Population")
if not rlayer.isValid():
    print("Layer failed to load!")
QgsProject.instance().addMapLayer(rlayer)

# Add Intensity
path_to_intensity = "..../mi.shp"
vlayer = QgsVectorLayer(path_to_intensity, "Intensity", "ogr")
if not vlayer.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(vlayer)

print(QgsProject.instance().mapLayers())
...............................................................

The base map is added as an .xml file, there is a population raster .tif file as well as an intensity shapefile. Ultimately, I run population statistics with the shapefile, to roughly work out the affected population within each intensity (as can be seen on the maps below). When I use the QGIS Python Console and the print composer, then exporting the layout as a .png file, the map looks like this and it all works perfectly.

However, when I use PyQGIS and add the data, run statistics, change symbology, create a new layout/print composer, add map features and then export the map, the exported map does not print successfully. I can see that the statistics and symbology has been completed but the raster base map is not seen on the map. I know that the map has been added to QGIS and the layout, because the blue water is displayed and if the file was not loaded properly, there would be white instead. I have tried to change the CRS of the layer but still no success. I have also tried moving the layers up/down but no success. The scale is wrong too but I think this is due to the base map not loading properly.
Does anyone know why the raster layer is not appearing in my map and why the export of the map does not include this layer using Python?
Could it potentially be down to the export resolution?
It is weird as I have another script that runs a similar process and has the same base map and it export perfectly.

A WeTransfer link to download the three files is - https://we.tl/t-xgrHBo8Fos


texts: wkid error using buffer GeometryService - ArcGIS JS API 3.2

texts:
First of all, going to this url for the geometry service referenced in the documentation shows it's been deactivated so I instead chose this url - both with /GeometryServer and /GeometryServer/buffer.

I suspect it has something to do with the spatial reference of my geometry or buffer param but I've used different GCS and PCS with the same result.

require(["esri/map",
     "esri/SpatialReference",
     "esri/geometry/Point",
     "esri/symbols/SimpleMarkerSymbol",
     "esri/symbols/SimpleLineSymbol",
     "esri/graphic",
     "esri/Color",
     "esri/tasks/BufferParameters",
     "esri/tasks/GeometryService"],
function(Map, SpatialReference, Point, SimpleMarkerSymbol, SimpleLineSymbol, Graphic, Color, BufferParameters, GeometryService) {
var map = new Map("map", {
    basemap: "streets",
    center: [-97.744483, 30.272705],
    zoom: 15
});
//Locate By Lat/Long
document.getElementById("btnLocate").addEventListener("click", function() {
    var latitude = Number(document.getElementById("latitude").value);
    var longitude = Number(document.getElementById("longitude").value);
    var pointGCS = new esri.geometry.Point(longitude, latitude, new SpatialReference({ wkid: 4326 }));
        var pointWMS = esri.geometry.geographicToWebMercator(pointGCS);
        var crossHair = new SimpleMarkerSymbol(SimpleMarkerSymbol.STYLE_CROSS, 12,
            new SimpleLineSymbol(SimpleLineSymbol.STYLE_SOLID, new Color([255, 0, 0]), 2), new Color([255, 0, 0]));
        var centerpoint = new Graphic(pointWMS, crossHair);
//Set buffer parameters
    var params = new BufferParameters();
    params.bufferSpatialReference = new SpatialReference({ wkid: 32662 });
    params.outSpatialReference = new SpatialReference({ wkid: 32662 });
    params.unit = GeometryService.UNIT_STATUTE_MILE;
    params.distances = [0.25, 0.50, 1];
    params.geometries = [centerpoint];
    var geoService = new GeometryService("http://tasks.arcgisonline.com/arcgis/rest/services/Geometry/GeometryServer");
    geoService.buffer(params);
    map.centerAndZoom(pointWMS, 15);
    });
});


texts: Make hangs when attempting to build zip code TopJSON file on Mac OS X

texts:
Thanks so much to all involved with the us-atlas project. I am learning how to use d3, leaflet and other tools to make maps and this is really helpful.

I have grabbed files from the us-atlas project on github (https://github.com/mbostock/us-atlas) and have had success doing the basic make for US Counties.

However, the make for US Zip codes seems to hang. I have tried running in background with & and also using nice to increase the priority but it still seems not to finish.

In response to Brad's question, the process has run in background for upwards of 10-11 hours overnight. The top command line utility at that point shows that node has 0% of CPU (whereas at the start I monitor and see it is taking 100% of CPU for the first hour).

Could it just be that the OS falls into sleep mode and it kills the processing?

Is there anything I can do to get OS X to keep the make (and subsequent node process) a priority and keep it going?

How long should this take with a 2.6GHz Intel Core i7 and 8MB RAM?

As an alternative I am mostly interested in MD Zip codes. Is there a way to restrict a TopoJSON build to just those zip codes?


texts: How to delete user defined project variables using Python

texts:
I am trying to delete a variable from QGIS Project variables. I didn't find the functionality in QgsExpressionContextUtils (I use it to set a variable).
I also tried to do this:

project_context_scope = QgsExpressionContextUtils.projectScope()
project_context_scope.removeVariable('test')

It does remove from project_context_scope. But, when I open the QGIS Project properties to see the list of variables, the deleted variable is still exist.

May be I was wrong in getting the instance of project_context_scope, since it's not a singleton and projectScope create new object.


texts: Finding an open data layer which includes towns , villages and ciities in the UK

texts:
I am looking for a data layer that i can import into Qgis to show towns, villages and cities. either of the three in an individual shp file would also be useful.

this is part of a project in which i am siting locations for windfarms and one critera for there posistioning is that they have to be over 2000m from a large dwelling

thus far i have found a buildings layer in Os vector district however this is not ideal.


texts: Where can download a DEM for whole world in a single low resolution image?

texts: I work in Blender to make Earth globe and I need elevation data for whole earth with low spatial resolution to make relief as a terrain.


texts: Applying graduated renderer to polygon using PyQGIS 3?

texts:
I'm trying to symbolize a polygon layer using a graduated color ramp, with no success.

In the attribute field of my layer I have a column: 'PatCNT' with integer values from 2 to 6 for each row.

I've been trying to adapt QGIS v2.x code from Applying graduated renderer in PyQGIS? or from QGIS python programming cookbook (also for QGIS 2.x)
I'm trying to use the examples from the PyQGIS cookbook describing how to use the graduated symbol renderer. I have pretty much copied the code from the cookbook example below and tried to adapt it to my case.

from qgis.PyQt import QtGui
from qgis.core import *

uri = '/Users/ep9k/Desktop/SandraMonson/TestZips.shp'

myVectorLayer = QgsVectorLayer(uri, 'testZipCodes', 'ogr')
myTargetField = 'PatCNT'
myRangeList = []
myOpacity = 1

#Make our first symbol and range
myMin = 2.0
myMax = 4.0
myLabel = 'Group 1'
myColor = QtGui.QColor(Qt.red)
mySymbol1 = QgsSymbol.defaultSymbol(myVectorLayer.geometryType())
mySymbol1.setColor(myColor)
mySymbol1.setOpacity(myOpacity)
myRange1 = QgsRendererRange(myMin, myMax, mySymbol1,
myVectorLayer.geometryType())
myRangeList.append(myRange1)

#now make another symbol and range
myMin = 4.1
myMax = 6.0
myLabel = 'Group 2'
myColor = QtGui.QColor(Qt.blue)
mySymbol2 = QgsSymbol.defaultSymbol(myVectorLayer.geometryType())
mySymbol2.setColor(myColor)
mysymbol2.setOpacity(myOpacity)
myRange2 = QgsRendererRange(myMin, myMax, mySymbol2, myLabel)
myRangeList.append(myRange2)

myRenderer = QgsGraduatedSymbolRenderer('', myRangeList)
myRenderer.setMode(QgsGraduatedSymbolRenderer.EqualInterval)
myRenderer.setClassAttribute(myTargetField)
myVectorLayer.setRenderer(myRenderer)
QgsProject.instance().addMapLayer(myVectorLayer)

I am just trying to have "group 1" in red and "group 2" in blue for now. Then I can customize it further in the future if I want.

Depending on what I change, I get one of many error messages...


texts:
I have a geoserver instance providing imagery, hillshade and OSM as WMS services. These tiles are mapped in Leaflet and the user can customize layer order, opacities etc.

I would like to be able to generate large format, high-res PDF from selected map layers and styles. I think this would be best completed server side.

Does anyone have suggestions for a server-side method for rendering pdf's from WMS layers?


texts: Using one world file for multiple rasters

texts: I have 400 PNG files all of the same area that are not gereferenced. I have gereferenced one of them and want to use the world file for all the others. I could copy it and rename the files but that will take forever. It there a way in QGIS to point the png files to a worldfile?


texts: Creating map with multiple selectable layers using ArcGIS

texts:
I want to make a map with an option for the user to choose between 3 layers to view. All layers with a certain type of point. If possible, I would also like the user to view more than one layer at a time if they select it.

Is this possible within ArcGIS?


texts: How to perform Intersection using ModelBuilder iteration?

texts:
I would like to take the intersection of two areas (one randomly generated), and do this n times. I can do it and get one result, but I can't save all the intersection areas to a file at the end.

The problem is that when I run my model more than once it says that my buffer output file already exists, which I can imagine would happen again for the intersection. What I'd like to happen is to run my model where I produce new random points everytime, then buffer them and intersect the dissolved buffer with another area and then get out the area. Then I want to do this several times.

If this isn't worded very well then just let me know and I can explain what I mean in another way. Thanks a lot


texts:
I want to get attributes from object what was selected on map. I understand how to get selected object:

Dim nrecords as Ineger
Dim obj as Object
nRecords = SelectionInfo(SEL_INFO_NROWS)
For i = 1 To nRecords
   Fetch rec i From selection
   obj = selection.obj

Now i want to get any attributes from this object. How can i do this?


texts: GRASS via Processing Toolbox doesn't work in QGIS 2.8.3 on Mac?

texts:
I am quite new to QGIS and working with the external applications such as GRASS.

When trying to launch a GRASS tool from the processing toolbox I get the following error:

  Missing dependency. This algorithm cannot be run :-(
  This algorithm
  requires GRASS to be run. Unfortunately, it seems that GRASS is not
  installed in your system, or it is not correctly configured to be used
  from QGIS
  Click here to know more about how to install and configure
  GRASS to be used with QGIS

When I am using GRASS from the Plugins menu it works. I have GRASS 6.4 installed to my system and running, GRASS 7 installed, but it doesn't run, probably due to some known El Capitan issue. Anyhow, GRASS 6.4 should be working. Moreover, there is also the GRASS that comes with the QGIS installation.

Unfortunately I messed up the path to GRASS that was originally set in Processing -> Options before I even tried if the original setup works. I tried really every possible path I could think of, to the GRASS-6.4 installation in my applications folder (/Applications/GRASS-6.4.app) and to the GRASS folder in the QGIS installation (/Applications/QGIS.app/Contents/MacOS/grass), I don't get the GRASS tools to working. I reinstalled QGIS to restore the original settings, but obviously I have to delete the particular settings before (but where?), because the path was still as I set it before reinstallation.

SAGA is working from the processing toolbox.

So, long story short:

Any suggestions what could be the issue? Anyone experiences with running GRASS via processing toolbox with QGIS 2.8 on El Capitan?
It would be already a great help if someone could give me the original path to GRASS set at "Processing" --> "Options"


texts: Difference between transformation values in Proj4 npm and epsg.io

texts:
The values are not matching with other online tools for the conversion using Proj string and in the case of wkt format, the conversion using -ive value doesn't match.
While using WKT string in epsg 25832 to convert to epsg 31467 there is a difference in the value when compared to epsg.io and this difference becomes more while moving towards the -ive boundary side. This difference is not seen while using the proj4 string.
EPSG 25832 to EPSG 31467 : X =-1877994.66 Y = 3932281.56
Result from WKT: 1121242.585 3933381.455
Result form proj4: 1121120.431 3933393.947
Result from epsg.io: 1121120.43 3933393.95


texts: ArcMap ignoring TiledLabelsToAnnotation in ArcPy script?

texts:
I'd write python script that execute TiledLabelsToAnnotation_cartography tool.
This script is call by a tbx file.

No problem when i launch my python script in Visual Studio Code, result is correct.
Now, i need to valid my script in ArcMap to share my geoprocessing as web service .sd file (arcgis server).
But arcmap ignore TiledLabelsToAnnotation_cartography tool and break my script... So, layers generated by this function missing and i have failed execution because i need this layers to realize a clip.

  /!\  Failed to execute. Parameters are
  not valid. ERROR 000732: Input Features: Dataset
  C:\Users...\scratch.gdb\LAYER_XXX_SEP150
  does not exist or is not supported

Already try :

realize TiledLabelsToAnnotation_cartography execution with python geoprocessing window with same parameters =  excecution success.
Manually import missing layers in gdb = execution success.

Why arcmap ignore this function ?

Code :

\# Entry : MXD_ETQ contain labels layers
MXD_ETQ = arcpy.mapping.MapDocument(path)

\#copy MXD to preserve original MXD from additionnal annotation layers
MXD_copy_path  = os.path.join(temp_folder, "_labels")
MXD_ETQ.saveACopy(MXD_copy_path)

\#tiled labels to annotation ---> Ignored by ArcMap
arcpy.TiledLabelsToAnnotation_cartography(
        MXD_copy_path, df_name , EXTENT_LAYER , TEMP_GDB , "ETQ", "SEP", 150,
        "", "", "", "", "STANDARD", "NOT_GENERATE_UNPLACED_ANNOTATION")


texts: Why do I get an error message at QGIS training manual 3.3?

texts:
I have a problem following the QGIS training manual lesson 3.3 (Ratio Classification).

You have to add a column to an SQlite layer. The new column is AREA. I followed the instructions and added it as shown, then saved it, and received the following error message (in orange on the screen)

  Commit errors: Could not commit changes to layer land use

Is the manual wrong, or am I doing something wrong?
Can anyone else duplicate what the manual asks?


texts: How to add a GeoWebCache gridset to all layers in GeoServer

texts:
I have 300+ layers in my GeoServer layer library, I have just created a custom gridset I would like to use for all of my layers. I don't want to go through my whole library adding the gridset manually through the UI, so I have played around with the GeoServer rest API and I am able to modify layers with the following command:

curl -v -u admin:geoserver -XPOST -H "Content-type: text/xml" -d @ui_UI00000ENNR.xml  "http://localhost:8080/geoserver/gwc/rest/layers/ui:UI00000ENNR.xml"

And this is the contents of the ui_UI00000ENNR.xml file:

<?xml version="1.0" encoding="UTF-8"?>
<GeoServerLayer>
  <id>LayerInfoImpl-288ccd00:163bb872e8d:-7ff6</id>
  <enabled>true</enabled>
  <inMemoryCached>true</inMemoryCached>
  <name>ui:UI00000ENNR</name>
  <mimeFormats>
    <string>image/png</string>
    <string>image/jpeg</string>
  </mimeFormats>
  <gridSubsets>
    <gridSubset>
      <gridSetName>EPSG:3857</gridSetName>
    </gridSubset>
    <gridSubset>
      <gridSetName>EPSG:4326</gridSetName>
    </gridSubset>
    <gridSubset>
      <gridSetName>EPSG:27700</gridSetName>
    </gridSubset>
    <gridSubset>
      <gridSetName>EPSG:900913</gridSetName>
    </gridSubset>
  </gridSubsets>
  <metaWidthHeight>
    <int>4</int>
    <int>4</int>
  </metaWidthHeight>
  <expireCache>0</expireCache>
  <expireClients>0</expireClients>
  <parameterFilters>
    <styleParameterFilter>
      <key>STYLES</key>
      <defaultValue></defaultValue>
    </styleParameterFilter>
  </parameterFilters>
  <gutter>0</gutter>
</GeoServerLayer>

My plan is to create a template xml file and then loop through my list of layers, populating the xml file and running the curl command for each layer. The only problem is I don't know the id for each layer, <id>LayerInfoImpl-288ccd00:163bb872e8d:-7ff6</id>, in the example above. Does anyone know how I could pull all layer names and there associated ids from GeoServer? Or is there a much quicker way of doing what I want to accomplish?


texts: Show QGIS feature count via console on existing layers

texts:
I want to make a script that toggles feature count on for all vector layers. I followed the following article: Show feature count of layer via Python console / PyQGIS.

Instead of adding new layers with feature count I'd like to switch on GUI feature count on existing layers. I adapted the code in the article so that myLayer is substituted with an existing layer that I retrieved by listing the layers.

The new layer made by the unadapted code to myLayerNode:

<qgis._core.QgsLayerTreeLayer object at 0x000000001FEC4950>

The existing layer after saving it in myLayerNode:

 <qgis._core.QgsLayerTreeLayer object at 0x000000001FEC49D8>

The new layer feature count can be toggled on/off using:

myLayerNode.setCustomProperty("showFeatureCount", True)

To my knowledge this should work but it doesn't...

#Make a list of layers
qgis.utils.iface.iface = iface
layers = qgis.utils.iface.iface.legendInterface().layers()

## reference to the layer tree
root = QgsProject.instance().layerTreeRoot()

## pick one layer to change featurecount
myLayer = layers[5]
myLayerNode = QgsLayerTreeLayer(myLayer)
#root.insertChildNode(5, myLayerNode)

## set custom property
myLayerNode.setCustomProperty("showFeatureCount", True)

myLayerNode.customProperty("showFeatureCount")
## the result is True

myLayerNode.customProperties()
## the result is a list [u'showFeatureCount']

Any idears?


texts: I have a property boundary and elevation data on a map I'm making. I would like to "trim" the elevation data so I can only see contours within the property boundaries not outside of them. How could I do this? See attached image.


texts: Does GeoServer's WorldImage plugin use UTM or WGS84 (EPSG:4326)?

texts:
I have some PNG images that I want to put on an OpenLayers map using GeoServer as a WMS, and it seems as though the easiest way to do that would be to create an accompanying .pnw file (ESRI world file) for each .png and load them using GeoServer's WorldImage plugin.

I also have the longitude, latitude coordinates of nine points on each rectangular image: {top, middle, bottom} x {left, middle, right}, and I can use five of them to convert pixel x,y into longitude, latitude x',y':

where (mm_lng, mm_lat) is the position of the middle-middle point, (mr_lng, mr_lat) is middle-right, (tm_lng, tm_lat) is top-middle, etc. for the middles of the four sides.  This is a linear approximation to the full coordinate transformation, minimizing error at the center of the image.

The World Image format describes an affine transformation in the same way: the six parenthesized terms above can be interpreted as the six lines of the .pnw file (in A, B, C, D, E, F order) if the map units are longitude, latitude coordinates (that is, WGS84 a.k.a. EPSG:4326).  The World file Wikipedia page says "World files describing a map on the Universal Transverse Mercator coordinate system (UTM) use these conventions..." and goes on to describe the UTM application, though it also says, "World files do not specify a coordinate system; this information is generally stored somewhere else."

If the result of the .pnw conversion is interpreted as UTM, then I would need to do additional work to transform from WGS84 into UTM, and GeoServer would also need more information to specify which of the 60 UTM zones that the .pnw converts into.  I don't see how to provide that information to GeoServer.  Am I right in assuming that GeoServer interprets the .pnw-converted coordinates as WGS84?


texts: Installing Data Interoperability when already authorized on ArcGIS Administrator?

texts:
If data interoperability is authorized on ArcGIS Administrator but is not yet installed, how does one do so in ArcGIS 10.2.2?

I have the license, now how do I get the extension? I have already checked the Customize menu and data interoperability is not listed. Nor can I find it doing a very thorough search in Windows Explorer. My license has not expired, but I cannot figure out how to install the darn thing!

I forgot to say that I am looking for the executable file that installs the extension. Online help has just said that it is in my "media file".


texts: Problem with algebra formulas for work with gdal_calc.py

texts:
I'm trying to update some values of a raster based on another raster.

1st raster (A): landuse data (values from 10 to 240)

2nd raster (B): mask data (only nodata values and a given value, for example: 20)

Objective: for a given pixel, if Pixel A == 212 and Pixel B == 20, replace pixel A with pixel B value.

So, my 1st command is: B*(A==212). That works fine, I'm getting a raster with pixels that matches the objectives and with the new value.

I'm having problems with the second raster manipulation. I dont know how to get a new full raster with values updated.

Anyone have an idea of how I can do that in, if possible, only one math formula ?


texts: Can I perform an analysis operation at a set time in ArcGIS Online?

texts:
I use ArcGIS Online to maintain data collected in the field for the company I work for. I perform an Overlay analysis operation at the end of each day when workers are done collecting. I perform the analysis manually each day. I was wondering if there is a way to set up the Overlay analysis function to automatically run each day at a specified time. That would prevent me from having to run it manually each day.

I could not find anything on the Internet saying this was possible. The closest thing I found was that I can use the "Rerun Analysis" tool on layers in ArcGIS Online. This would still involve me manually performing the analysis, but would speed things up a little possibly.


texts:
I need to join two layers (a fishnet grid and land cover). I'd like the join to occur only if the the land cover occupies 50% or more of each fishnet grid (See image below).

Can anybody provide an example or point me in the right direction?

More details:
Ideally I only need the grids to be created that contain 50% or more of the land cover  layer.


texts: Topographic Index - Why is upslope contributing area divided by downslope contour width and not upslope contour width?

texts:
In the Topographic Index calculation, a is the upslope contributing area per unit contour length.  Is the "unit contour length" calculated from examining the contour widths of the downslope cells that water in the current cell flows to?  Or is it the contour widths of the upslope cells that flow to the current cell?

In Quinn and Beven's work, it seems to be the the downslope cells.  Why?  How does that make more hydrological sense than using the upslope cells' contour widths?  Can some justify the use of the downslope (from the current cell) cells' contour widths vs the upslope (from the current cell) cells' contour widths for me?

Edit - Hopefully these pictures and this description provide the details people are looking for.

Let's examine a 3 x 3 grid with elevation values and a 10m cell size resolution.  Let's assume we are working with a MFD like algorithm.  The grid is pictured below:

In this grid, the center cell, with a value of 5, receives flow (and therefore area) from the upslope cells with higher elevation values.  So it receives flow (and therefore area) from the cells that contain elevation values of 8, 7, and 6.  Different algorithms (Quinn, Freeman, etc.) can dictate how much area the upslope cells pass to the center cell based on the fraction that particular slope (upslope cell to center cell) contributes to the sum of all the slopes (upslope cell to all downslope cells).  Anyway, however it is decided, accumulated area ends up in the center cell passed to it from its upslope neighbors.  This accumulated area can also be called the upslope contributing area.

In this image we see the downslope cells that the center cell passes area/flows into.  Now in Quinn's description of his algorithm, he says :

  Thus the actual value of a/tan ß for the current cell will be the
  total area draining through that cell acting across a contour length
  equal to the sum of the downhill contour lengths and divided by the
  weighted average tan ß value.....(THE PREDICTION OF HILLSLOPE FLOW PATHS FOR DISTRIBUTED HYDROLOGICAL MODELLING USING DIGITAL TERRAIN MODELS)

What I am confused about is this:

Since the upslope contributing area is dictated by the both the slope between the upslope cells and the center cell as well as contour width determined by the direction (cardinal vs. diagonal), why are the downslope contour widths used?  The upslope contributing area per unit contour width, the numerator in the TWI/TI, seems to dictate how wet an area will be due to accumulated flow from upslope cells - why divide it by the the contour widths of the flows being passed downslope?  And the denominator seems to get at how wet a cell will be in relation to it's steepness, with the steeper the drop to the downslope cells the larger the sum of their slopes, and the larger the denominator will be and the lower the overall TI value will be.

The denominator I understand fine, but I can't seem to wrap my brain around why the area passed into the cell is divided by the downslope contour widths....and not the contour widths of the flows entering the center cell from the upslope cells.


texts: Can't use MapServer to create custom basemap or ArcGISTiledMapServiceLayer

texts:
I have a NZ MapServer, that has "Single Fused Map Cache: true", yet I am unable to use it for custom basemap or ArcGISTiledMapServiceLayer.

This works:

esriBasemaps.customBasemap = {
  baseMapLayers: [{url: "https://services.arcgisonline.com/ArcGIS/rest/services/Specialty/DeLorme_World_Base_Map/MapServer"}],
  title: "customBasemap"
};

var map = new Map("map", {
  basemap: "customBasemap",
  center: [172.62982, -43.523158],
  zoom: 17
});

This doesn't:

esriBasemaps.customBasemap = {
  baseMapLayers: [{url: "https://services.arcgisonline.co.nz/arcgis/rest/services/Imagery/newzealand/MapServer"}],
  title: "customBasemap"
};

var map = new Map("map", {
  basemap: "customBasemap",
  center: [172.62982, -43.523158],
  zoom: 17
});

I've seen a couple of people suggesting to use ArcGISTiledMapServiceLayer instead, but that would mean I draw a map and then another on top of it, so I'd like to use it as a Basemap.

Can someone please tell me what I am missing?


texts: Converting a shp file to a list of latitudes and longitudes

texts: I have a shp file with the latitudes and longitudes for the boundaries for various geographic regions.  I'd like to convert the shp file into lists of just latitudes and longitudes for the regions.  I'll then use these lists in my program to determine which region a specific point resides in.  Is there a quick way to do this conversion or will I have to read up on the APIs?  Even converting the shp file into GML would be a help as I could easily access that programmatically.


texts: How to link a theme to a map layout in QGIS?

texts: I have different map layouts In QGIS which use different map themes. Currently, I have to choose a theme first before opening its corresponding map layout. Is there a way to link a theme to a map layout, so that whenever I open a map layout, its theme will be turned on and the map layout is showing what are desired automatically?