texts: Leaflet distinguishing sublayers within one GeoJSON layer

texts:
I would like to have more than one layer within one existing GeoJSON layer.

I found some example here:

http://www.gistechsolutions.com/leaflet/DEMO/sports/sports.html

which refers to one .json file and afterwards the features are collected from the properties.

I was trying on my code do something more simple, as I want to have the same colour, but only smaller radius of my icon.

My GeoJSON code looks like this:

 var sitec = {
    "type": "FeatureCollection",
        features: [
            {
             type: "Feature",
             "properties": {
                 "Title": "Sitec IS",
                 "Head": "7400 Beach Drive",
                 "Description": "Gavin Sinclair",
                 "Value":2
        },
            "geometry": {
            "type": "Point",
            "coordinates": [
            0.16964435577392578,
            52.29220753602784
                ]
                  }
                    },
 {
  "type": "Feature",
  "properties": {
"Title": "JHG",
    "Head": "Shortstanton sidings",
"Description": "Conor Murphy",
"Value":1,
"URL": "JHG",
"Pict": "image.png"
 },
  "geometry": {
    "type": "Point",
    "coordinates": [
      0.05458831787109375,
      52.29163006501503
    ]
  }
 }
              ]
        };

and what I would like to achieve - set the file icon size based on the propeerties.Value   but keep it in the same group.

My marker customisation looks like this:

var geojsonMarkerOptions = {
radius: 8,
fillColor: "#ff7800",
color: "#000",
weight: 1,
opacity: 1,
fillOpacity: 0.8
};

and again, I would like to have everything the same, apart from radius, so just in case I set:

var geojsonMarkerOptions3 = {
radius: 2,
fillColor: "#ff7800",
color: "#000",
weight: 1,
opacity: 1,
fillOpacity: 0.8
};

My code looks as follows:

 var sitis = L.geoJSON(sitec, {
pointToLayer: function (feature, latlng) {
feature.properties.myKey = feature.properties.Title + ', ' +
 feature.properties.Head
feature.properties.Value
label = String(feature.properties.Title)
    return L.circleMarker(latlng, geojsonMarkerOptions).bindTooltip(label,
  {permanent: true, direction: "center", className: "my-
  labels"}).openTooltip();
  return L.circleMarker(latlng, geojsonMarkerOptions3);
},
    onEachFeature: function (feature, layer) {
        layer.bindPopup("
"+feature.properties.Title+"

Address: "+feature.properties.Head+"

"+feature.properties.Description+"

 Website:"+feature.properties.URL+"

");

  }
    })
    .addTo(map);

Map remains the same and console says nothing.

Does anyone knows how to do it? Or shall I write all this stuff separately, creating sth like completely new layer?


texts: Problem with multi-threading and geoprocessor tool in c#

texts:
What I am doing: I am trying to run a geoprocessor tool in a background thread other than UI thread. I am passing string path of arcObjects and creating that arcObject in my background thread. Then I am passing those arcObjects to geoProcessor tool.

In UI thread, Button's on click method:

private void btnFindInterction_Click(object sender, EventArgs e)
{
    List<IFeatureLayer> layers = GetAllMapFeatureLayers();

    string featureLayerPath1 = GetFeatureLayerPath(layers[0]);
    string featureLayerPath2 = GetFeatureLayerPath(layers[1]);

    Thread worker = new Thread(new ThreadStart(() =>
        {
            IFeatureClass featureClass1 = GetFeatureClassFromPath(featureLayerPath1);
            IFeatureClass featureClass2 = GetFeatureClassFromPath(featureLayerPath2);

            var result = RunIntersectTool(featureClass1, featureClass2);
            OnCalculationComplete(intersectedFeatureClasssPath);
        }
    ));

    worker.SetApartmentState(ApartmentState.STA);
    worker.Start();
}

The RunIntersectTool method:

string RunIntersectTool(IFeatureClass featureClass1, IFeatureClass featureClass2)
{
    IWorkspace workspace = GetTempGDBWorkspace();
    string intersectFeatureClassName = "Intersect_result";
    string intersectFeatureClassPath = workspace.PathName + "\\" + intersectFeatureClassName;
    try
    {
        var gp = new Geoprocessor { OverwriteOutput = true, AddOutputsToMap = false };
        IGpValueTableObject gpValueTableObject = new GpValueTableObjectClass();
        gpValueTableObject.SetColumns(1);
        var inFeature1 = imperviousFeatureClass as object;
        gpValueTableObject.AddRow(inFeature1);
        var inFeature2 = watershedFeatureClass as object;
        gpValueTableObject.AddRow(inFeature2);

        var intersectTool = new Intersect
        {
            in_features = gpValueTableObject,
            out_feature_class = intersectFeatureClassPath
        };
        gp.Execute(intersectTool, null);
        return intersectFeatureClassPath;
    }
    catch (Exception e)
    {
        return null;
    }
}

What problem am I facing:
The tool runs perfectly in some machines, but gives error in some machine: "000732 : : Dataset  does not exist or is not supported" And this error never happen if I don't use threading. I read help files of Running a geoprocessing tool using background geoprocessing, but I couldn't find anything that says I can't pass arcObjects.

There is another problem I am facing. Whenever I ran the tool, it works fine first time(in some machines), but gives error when it creates an instance of geoprocessor at second time. Error message: "Attempted to read or write protected memory". I guess this error is also related to multithreading and arcObjects.


texts: overlay my shapefiles stored in database to the web/browser

texts: I am very much strrugling on how to display my shapefiles stored in postgis database to the browser or to the web.. is there anyone there knows how to do it.? I am using geodjango framework in python,. I already achieve the displaying of maps in admin interface but when I create my public page for displaying the maps, its make my head aching,.I did not konow how to overlay my shapefiles stored in database to the base map openlayers / openstreetmap.


texts:
I am trying to show ECW ortophoto as base layer in openlayers.
So i found this thread on GIS SE:

    Geoserver 2.1.1 with ECW (v3.3) support

1. Download Geoserver 2.1.1 from http://downloads.sourceforge.net/geoserver/geoserver-2.1.1-bin.zip 1a. extract to c:\bin\geoserver

2. Download imageio 1.1.1 installer from http://java.net/projects/imageio-ext/downloads/download/Releases/ImageIO-Ext/1.1.x/1.1.1/windowsInstaller/windows32-imageio-ext-installer-gdal-mrsid-ecw-1.1.1.zip 2b. Install it to c:\imageio 2c. Remove c:\bin\geoserver\webapps\geoserver\WEB-INF\lib\imageio-ext-1.0.8 2d. Copy c:\imageio\lib*.* to c:\bin\geoserver\webapps\geoserver\WEB-INF\lib (imageio installer will also copy gdal dll, ecw sdk dll to your java\bin directory)

3.Add Environment variables for GDAL_DATA c:\imageio\gdaldata GDAL_DRIVER_PATH c:\imageio\gdalplugin

4.Start geoserver

I did all this steps except number 3, what its means to Add Environment variables?
How I do that?


texts: Install the 5.11.2 version of PyQt5 using OSGeo4W setup

texts:
How can I install the 5.11.2 version of PyQt5 using the OSGeo4W setup in Windows10? The default version of PyQt5 installed from OSGeo4W is 5.11.3 which is different from the QGIS compilation version, as in the pictures below.

My problem is when I run pylupdate5 translation.pro to create a .ts file for internationalization of code I give this message:
Traceback (most recent call last):
  File "C:\OSGeo4W64\apps\Python37\lib\runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "C:\OSGeo4W64\apps\Python37\lib\runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "C:\OSGeo4W64\apps\Python37\lib\site-packages\PyQt5\pylupdate_main.py", line 21, in <module>
    from PyQt5.QtCore import (PYQT_VERSION_STR, QDir, QFile, QFileInfo, QIODevice,
ImportError: DLL load failed: The specified module could not be found.

So, after that, I manually uninstall PyQt5 5.11.3 and install the 5.11.2 version and run pylupdate5 translation.pro again and it works, but when I start QGIS this message appears:
Error message on launch:
Couldn't load SIP module.
Python support will be disabled.

How can I solve this?


texts: How to reproject Network Dataset from GCS_WGS_1984 to Web Mercator?

texts:
I need to know how to change the coordinate system of my feature dataset with its network dataset from GCS_WGS_1984 to Web Mercator.

Once I know how to do this, I would like to make a tool in Python or ModelBuilder because I need to do it many times.


texts: Why does the Esri-resource-proxy log users out of Esri Geoportal Server?

texts:
I'm seeing a problem whereby passing a query to the Esri Geoportal Server's REST API, via the Esri Resource Proxy, is causing users to be logged out of the Geoportal.

I can reproduce the problem like this:

install the Esri Geoportal Server on Windows using IIS
install the Esri Resource Proxy on the same machine
log in to the Geoportal (eg http://localhost:8080/geoportal/catalog/identity/login.page)
in a new tab in the same browser, send a query to the Geoportal REST API via the Esri proxy (eg http://localhost/PHP/proxy.php?http://localhost:8080/geoportal/rest/find/document?f=pjson)
switch back to the Geoportal tab, refresh the page, and note that the user is now logged out

This error does not occur for me using the sandbox Geoportal at http://gptogc.esri.com/geoportal/catalog/main/home.page (you can test it with this query)

Does anyone know why this is occurring? Is the problem occurring because the proxy is on the same machine as the Geoportal Server? Are there any known workarounds?

I don't know enough about how the login session works in the Geoportal Server, or what the Esri Resource Proxy is doing, to know how to debug this.


texts:
I have a script that buffers a selected point, and intersects that buffer with features that the user inputs. I would like to convert the table outputted from the Intersect process into a excel file. This is what I have so far:

import arcpy

#Set to current mxd and dataframe
mxd = arcpy.mapping.MapDocument ('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]

# Set overwrite option
arcpy.env.overwriteOutput = True
arcpy.env.workspace ="C:\Users\pierrej\Desktop\GIS Data"
wrk=arcpy.env.workspace
arcpy.AddMessage("Workspace added")

try:
# Get the input parameters for the Selection Tool
FClass = arcpy.GetParameterAsText(0)
Field = arcpy.GetParameterAsText(1)
Feature = arcpy.GetParameterAsText(2)
outPath = arcpy.GetParameterAsText(3)
bufferDistance = arcpy.GetParameterAsText(4)
IntersectFeature = arcpy.GetParameterAsText(5)
SpecifiedSites = arcpy.GetParameterAsText(6)

where_clause = """{} = {}""".format(arcpy.AddFieldDelimiters(FClass, Field),Feature)   # I have added extra pair of single quotes for selecting strings
arcpy.AddMessage(where_clause)

# Select the site
arcpy.SelectLayerByAttribute_management(FClass,"NEW_SELECTION", where_clause)

# write selected features to a new featureclass
arcpy.CopyFeatures_management(FClass, "SelectionSites")

# Report a success message
arcpy.AddMessage("Site selected!")

#Zooming to a selection set for the specific layer
df.zoomToSelectedFeatures()
df.scale = 300000
arcpy.RefreshActiveView()

# Run the Buffer tool
AreaOfInterest = arcpy.Buffer_analysis(FClass, outPath, bufferDistance)

# Report a success message
arcpy.AddMessage("Buffer made!")

#Clear Selected features
arcpy.SelectLayerByAttribute_management(FClass, "CLEAR_SELECTION")

# Report a success message
arcpy.AddMessage("Selection cleared!")

# Run the Intersect analysis
List = arcpy.Intersect_analysis([IntersectFeature, AreaOfInterest], SpecifiedSites, "ALL", "", "INPUT")

# Report a success message
arcpy.AddMessage("Intersect analysis made!")

#Process: Export to CSV table
arcpy.TableToExcel_conversion(List, wrk, "Table.csv", "COMMA")

# Report a success message
arcpy.AddMessage("Results exported to a Excel file!")

except:
 # Report an error messages
 arcpy.AddError("Could not complete")

 # Report any error messages that the tool might have generated
 arcpy.AddMessage(arcpy.GetMessages())

So I'm trying a couple of things and I keep getting this error message.
ERROR 000732: Output Location: Dataset C:\Users\pierrej\Desktop\GIS Data does not exist or is not supported. Even if I change the output path to a folder without spaces I get the same thing.


texts: Getting Layer DataSource information in ArcObjects with C#?

texts:
I have a mxd file with 6 layers. I want to get datasource information for each layer. Please see the code below. I want to know how exactly i should get layer's datasource information

IDocument doc = ArcMap.Application.Document;
IMxDocument mxDoc = doc as IMxDocument;

IMap map = mxDoc.FocusMap;

ILayer layer = map.get_Layer(0);// Layer[0];


texts: Appending data in PostGIS database with QGIS?

texts:
I 'm running postgis with qgis 2.6. I created a database and tables that i need to update with new datas provided with shapes. I cannot find how to insert them in table. All attributes both shapes and table are same. I don't despair to have an "easy" way instead of using an insert to command that i'm not sure to use correctly in postgis command.
I've tried the database manager but i only succeed in erasing all previous data in the table >_<

Window 7 64 bits
Qgis2.6
Postgis1.5
Postgrsql 9.4


texts: Maximum number of layers in OpenLayers2 with GeoExt

texts:
I have been developing web gis application that requires around 45 different layers served from GeoServer. I had created OpenLayers WMS layers and added them to GeoExt, everything works fine till the 35th layer after I add 36th layer then OpenLayersgives pink tiles.
Here is my snippet:

 WGS84 = new OpenLayers.Projection("EPSG:4326");
             WGS84_google_mercator = new OpenLayers.Projection("EPSG:900913");

             map = new OpenLayers.Map({
                 height: '300px',
                 // allOverlays:true,
                 width: '500px',
                 projection: WGS84_google_mercator, displayProjection: WGS84,
                 maxExtent: new OpenLayers.Bounds(80.0509262084961, 26.348379135131836, 88.20467376708984, 30.4458065032959).transform(new OpenLayers.Projection("EPSG:4326"), WGS84_google_mercator),
                 //  )
                 controls: [
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.ArgParser(),
        new OpenLayers.Control.Attribution()
                 ]
             });
nnt_dyear_child_layer=new OpenLayers.Layer.WMS("Child Death Rate Yearly",
                         "hostname1/geoserver/wms", {
                             name: 'nnt_indicator_d',
                             layers: "vpdms:nnt_dyear-child", transparent: true,
                             //viewparams: a,
                         }, {
                             layerid: 'nnt_dyear-child_layer',
                             singleTile: true,
                             visibility: false,
                             group: "Indicators/NNT"
                         });
 nnt_dweek_child_layer=new OpenLayers.Layer.WMS("Child Death Rate Weekly",
                         "hostname1/geoserver/wms", {
                             name: 'nnt_indicator_d',
                             layers: "vpdms:nnt_dweek-child", transparent: true,
                             //viewparams: a,
                         }, {
                             layerid: 'nnt_dweek-child_layer',
                             singleTile: true,
                             visibility: false,
                             group: "Indicators/NNT"
                         });
map.addLayers([nnt_dweek_child_layer,nnt_dyear_child_layer]);

and finally  a geoext tree to hold the layer list:
tree1 = new GeoExt.ux.tree.LayerTreeBuilder({
                 region: "west",
                // title: "Layers",
                 renderTo: 'layers',
                 width: 240,
                 height:200,
                 autoScroll: true,
                 rootVisible: false,
                 lines: false,
                 enableDD :true,
                 // widget custom properties
                 wmsLegendNodes: true,
                 vectorLegendNodes: true
             });

I have used google base layers. I'm using openlayers 2.11 and geoext 2
Btw I have already read Here and Official Here
Any help would be great.


texts: How to rename an attachment when using ArcGIS Collector for Android?

texts: I am using ArcGIS Collector for Android, v 18.0.3 Build 1033 and I don't see a way to rename attachments in the field.  Is this possible from within the Collector application?  Or do I need to save photos to the device and rename them, and then attachment them to the feature in ArcGIS Collector?


texts: Basemaps not available in my ArcGIS 10.6.1

texts: My ArcGIS 10.6.1 has only 6 basemaps, it doesn't have the gray basemap etc.


texts: In what coordinate systems do true circles appear round?

texts:
I have been developing in Google Maps for awhile and am switching to OpenLayers but am stuck on one point: I am drawing vectors that are mostly circular but also involve some custom geometry so I can't use createRegularPolygon. I'm creating an array of points and sending it to a OpenLayers.Layer.Vector object. This works great but the resulting 'circles' are only round near the equator, not at higher latitudes. They are squashed about the same amount as the basemap is at higher latitudes, which seems to makes sense. If I use Google as a base layer then they are the correct shape at all latitudes, but I'm making an offline mobile solution so that won't do.

What do I use for a base layer (or what parameters to set) so that my geometry remains round at all latitudes without using Google or online maps?


texts: Start session edit on a specific layer on the map?

texts:
I created a GIS using ArcObjects, my application needs to be connected with a geodatabase stored in SQL Server. I used ARCSDE to connect my application to the database and when the database is connected, my application loads two tables - the first one is called Commune and the second one is called Parcelle. Their esriDatasetType is:
esriDatasetType.esriDTFeatureClass.

A Parcelle belongs to only one Commune; and a Commmune can contain many Parcelles.

What I want to do is start the session and edit the Parcelle layer without the possibility to edit the Commune layer. When I use my code, I can edit both layers in the map:

public static void Start(IEngineEditor _engineEditor, IMapControl4 map, string targetLayerName)
{
    //If an edit session has already been started, exit.
    if (_engineEditor.EditState == esriEngineEditState.esriEngineStateEditing) return;

    //Start editing the workspace of the first feature layer found.
    for (int layerCounter = 0; layerCounter <= map.LayerCount - 1; layerCounter++)
    {
        ILayer currentLayer = map.get_Layer(layerCounter);
        if (currentLayer.Name != targetLayerName) continue;
        try
        {
            if (currentLayer is IFeatureLayer)
            {
                IFeatureLayer featureLayer = currentLayer as IFeatureLayer;
                IDataset dataset = featureLayer.FeatureClass as IDataset;
                IWorkspace workspace = dataset.Workspace;
                _engineEditor.StartEditing(workspace, map.ActiveView.FocusMap);

                ((IEngineEditLayers)_engineEditor).SetTargetLayer(featureLayer, 0);
                break;
            }
        }
        catch (Exception ex)
        {
            Dialog.info("Info", ex.Message);
        }
    }
}

My function has 3 parameters: IEngineEditor, IMapControl4 and the name of the layer that I want to edit, which is the layer Parcelle in my case.


texts: GeoKettle and scheduling transformations / jobs?

texts:
Is there any possibility to schedule a transformation or job with GeoKettle?

I would like to perform the same task every day.


texts: Checking where values are the same between 2 rasters

texts:
I am working in ArcGIS 10.3 and I am trying to discover where the output of two different models agree.

two models were run for the risk of contamination with a carcinogenic toxin. The output of each model was a raster with risk scores of 0, 1, 2, and 4 as the values of each grid cell. the two rasters show some areas with the same score and sometimes the two rasters have different scores for the same area. Is there any way in ArcGIS to analyse the rasters to determine which areas have the same scores and which areas that have different scores on both maps? Im trying to determine areas of uncertainty between the models outputs.


texts: Getting 10 days period, max NDSI from MODIS/Terra 500m resolution daily product(MOD09GA)?

texts:
I want to make a NDSI with a ten days period from daily MOD09GA with Google Earth Engine.

The code is a adaptation from this code from Nicholas Clinton at Maximum ndvi for 3 months average in Google Earth Engine

When I use the code (the comment in the code) for renaming images after their dates, all goes well but when i try to download it, gives me  error. The error is

  Error: Expected a homogeneous image collection, but an image with
  incompatible bands was encountered: First image type: 1 bands
  ([2001_01_21]). Current image type: 1 bands
  ([2001_01_22]). Image ID: 2001_01_22 Some bands might require explicit
  casts.

How can I rename those NDSI's with the date name so that I can export it right?

var pol3=ee.FeatureCollection('users/georgeboldeanu20/front_WGS');
var batch = require('users/fitoprincipe/geetools:batch')

var NDSI = ee.ImageCollection('MODIS/006/MOD09GA').map(function(image) {
  return image.select().addBands(image.normalizedDifference(['sur_refl_b04','sur_refl_b07'])).rename('NDSI').clip(pol3);
});
//This is the rename code
////var NDSI_name = NDSI.map( function(img){
/// var id  = img.id()
 ///img = img.select('NDSI').rename(id);
///return img.copyProperties(img,['system:time_start','system:time_end']); });
//print(NDSI_name);

var COLlection = function(collection, start, count, interval, units) {
  var sequence = ee.List.sequence(0, ee.Number(count).subtract(1));

var originalStartDate = ee.Date(start);

  return ee.ImageCollection(sequence.map(function(i) {
    var startDate = originalStartDate.advance(ee.Number(interval).multiply(i), units);

    var endDate = originalStartDate.advance(
      ee.Number(interval).multiply(ee.Number(i).add(1)), units);

    return collection.filterDate(startDate, endDate)
        .reduce(ee.Reducer.max());
  }));
};

var days= COLlection(NDSI_name, ee.Date('2001-01-01'),10, 10, 'day');

var test = ee.Image(days.first());
print(test)
print(days)

batch.Download.ImageCollection.toDrive(luni, "Folder", {region:pol3.geometry().bounds(),crs:'EPSG:4326', scale:300});


texts: Does anyone know if it's possible to run a moving window analysis on either a raster, vector, matrix or data frame of spatial data to calculate Fragstats class statistics?


texts: Notifying user before overwrite by Copy tool of target with same name in ArcGIS ModelBuilder?

texts:
At ArcGIS Desktop (advanced) 10.3.1, I've created a model that simply copies a file geodatabase.  Both the source and target locations are parameterized.  Here's a screen shots of the model:

And here it is in action:

Here's my problem: the tool overwrites any pre-existing .gdbs with the same name in the target path.  I want to modify the tool so that it notifies the user prior to overwriting.  How can this be done?  Nothing in the Copy tool instructions or environment settings solves the problem.


texts: Update shapefile attributes after update of data source

texts: In a project the data is stored in a relational database (so it could be moved to PostgreSQL) for presentation in a web repository to be presented in diagrams, tables, etc. During the last months columns with spatial information were added to some entities and the idea of presenting these over a map came up. GeoServer and PostGIS are now set up (with the PostGIS database being seperate from the repositories DB on another server), delivering a shapefile to the clients. There will be updates to the dataset of the web app about once a week. Therefore the .dbf file of the shapefile needs an update too. Now, I could write a trigger which in its final consequence creates a new/updated shapefile - but: What is an established way of handling & automating this case, where data is stored in one DB and the shapefile is held in another?


texts: Why is gdalwarp flipping pixel size sign?

texts:
I am using the following command:

gdalwarp -co COMPRESS=LZW \
    -ts 1024 1024 \
    -cutline cutline.json \
    -crop_to_cutline \
    input.tif \
    output.tif

the output raster has the same number of pixels as the input one, but is just cropped to a cutline.

The input raster has pixel size (from gdalinfo):

Pixel Size = (159.847049179404479,153.365855104948253)

but the output raster has:

Pixel Size = (159.820436280028900,-153.866089899794588)

Why is gdalwarp flipping the pixel size?

What can I do to have positive size everywhere?

This is with GDAL 2.1.2


texts: PostGIS: Query using spatial join using pgadmin3

texts:
I have created three tables ch01.restaurants (SRID 2163;columnsfranchise, geom), ch01.lu_franchises (lookup table, no geom field;columnsid, franchise), ch01.highways (SRID 2163;columnsgid, feature, name, state, geom) in PostGIS that I want to use for a query which includes a spatial join. Running the query in pgadmin3, the following happens:

Total query runtime: 32 msec
0 Zeilen geholt.

My problem is that the table that will be created doesn't contain any data. I get a table with the heading

franchise | total

Here's my code:

SELECT f.franchise, COUNT(DISTINCT r.id) As total
FROM ch01.restaurants As r
INNER JOIN ch01.lu_franchises As f ON r.franchise = f.id
INNER JOIN ch01.highways As h ON ST_DWITHIN(r.geom,h.geom,1609)
GROUP BY f.franchise
ORDER BY total DESC;

How can I alter this code so that my result will contain all the franchise restaurants within one mile of a highway?


texts: How to convert unix time to qt time type?

texts: How can I convert unix timestamp to qt time type using qgis field calculator? I would like to avoid importing my table into excel and convert the timestamp to qt time type.


texts:
I am doing accessibility analysis. I created network dataset in ArcMap 10.6
I made a fishnet and I typed 100 and 100 for cell width and height, leaving rows and columns blank. I continue with OD Cost Matrix. I can create origins and destinations. When I click solve in network analyst, after 20 minutes or more, solver fails. This happens right after writing lines. Some of you might say it could be RAM problem. My computer's RAM is 16 GB. Interesting enough, when I choose cell width and height as 500 in fishnet, I can't see any problem, it works fine. But I need to create cell width and height of the fishnet as 100. Size of the area is 604161642.987255 (shape_area).

What should I do?

I am new to network analysis therefore accessibility.


texts: Avoiding pop-up when having multiple hyperlink layers in ArcGIS Desktop?

texts:
In ArcMap 10.0, I have a layer of country polygons and a layer of city points. Both layers have a hyperlink setup, and it all works as expected. The only drawback is that when I select a city, it comes up with a pop-up asking whether I want to go the hyperlink of that city or to the hyperlink of the underlying country.

Is there any way of avoiding this pop-up, so that it goes to the city hyperlink when I click within x pixels of a city, and goes to the country hyperlink when I click elsewhere within the polygon.

I recently upgraded to ArcGIS v10.2, but the same functionality still exists. And based on the responses so far, I probably didn't explain my question very well.

First of all, I am not talking about a bug or error. The functionality is as expected, but I am looking for a way to overrule the default functionality.

Secondly, for simplicity I used countries and cities, but in reality it was contaminated sites (polygons) and borehole locations (points). The problem however doesn't change.

Below is an ArcGIS screenshot, showing the hyperlinked points inside the hyperlinked polygons.

When I click inside the polygon away from a point, it goes straight to the polygon's hyperlink. This is what I want.
When I click on a point that is not inside a polygon, it goes straight to the point's hyperlink. Again, this is what I want.
When I click on a point inside the polygon, it comes up with the hyperlinks selection box shown in the screenshot, asking whether I want the point's or the polygon's hyperlink. This what I expect and it is conform the specifications, but its NOT what I really want. I want it to go straight to the point's hyperlink and ignore the underlying polygon, because the point layer is the uppermost layer of the two.

Is this easily doable, or would I have to design a complex Python function? Given that this is such a negligible inconvenience, it's not worth to alter this. But if there's a quick tip, I'd happy to hear.


texts: postGIS st_equals function not working as expected

texts:
I have a situation where I need to find out if a geometry being imported into a postGIS 2.1 database matches one drawn previously (i.e we want to know if the imported geometry is less up-to-date than the existing one).

I thought I could use st_equals for this, but it isn't working as expected.

This is the code I am using to find if the new geometry exists:

SELECT id FROM "GeomArchive" WHERE ST_Equals(geom, ST_GeometryFromText(wkt,4326));

however, it always returns false even when I pass in a geometry extracted from the "GeomArchive" table.

I saw this question, but the solution didn't help.

I find "=" works, which seems to make it unlikely that the problem is with rounding, since I would have expected a false result with this as well.

I'm not sure what "=" compares. Something I read in the pg documentation seemed to indicate that it actually only compares bounding boxes.

Why is st_equals not working in this situation?


texts: Export/Convert Layer to Shapefile within ModelBuilder

texts:
I am finding that more and more I am being sent KML files to be converted into shapefiles to be used for mapping/analysis etc. I am finding it monotonous using the Export KML to Layer tool, then re-exporting the layer (right clicking on layer in ToC and exporting as shapefile.) This seems like a lot of unnecessary work and creates a lot of the dreaded 'temp' files.

I would like to create a model to speed up the process. I am getting stuck however finding a way to export the output layer (from running the Export KML to Layer tool) to a shapefile.

Is there a way to do this? Or alternatively, export the .gdb the is also exported from running the Export KML to Layer tool?

I am also open to suggestions on how to do this using Python!


texts: Using WATERSHED to delineate tree crowns in ArcGIS Desktop?

texts:
I'm trying to delineate and count trees from a DSM/DTM generated from a UAV survey. It's from a citrus grove where trees are growing in different spacing. I have developed a workflow that works for most trees, but has some false positives. I'm stuck with reducing these false positives.

I'm using the concept of detecting local maxima and subsequently treat the trees as watersheds to delineate tree crowns. I'm using ArcGIS, for training purposes I'd like to only use this.

I'll provide my workflow below, additionally you can see some images with comments here: http://imgur.com/a/fLjia

Workflow:

extract tree height by subtracting the DSM from the DTM
clean up the resulting tree height raster with CON to remove negative values
use low pass filter to smooth the raster a bit (FILTER)
invert the raster. Trees are now holes. (RASTER CALCULATOR)
FOCAL FLOW to detect local minima (Value = 255)
Extract local minimas, use CON to replace 255 with the actual tree height from inverted tree height raster. I now have a raster which only has pixel values at the tree tops and NoData inbetween.
Convert tree tops to point features (= Input for pour point in WATERSHED)
Optional: Remove identified trees below a certain height (I used 1 m as cut-off)
IsNull to generate a raster with 1 and 0, where tree tops are 0 and the space inbetween is 1.
SetNull to use this 0-1 raster to replace the pixel-values of the tree tops in the tree height raster with NoData (this will be the pour point for WATERSHED)
Use FILL to remove sinkholes
FLOW DIRECTION to generate the flow direction
WATERSHED to generate the water sheds (= the tree crowns)
Invert the raster back to normal
Remove areas around trees (grass, soil, etc) by removing all values below a certain height (I used 1 m again).
This workflow works fine for most trees.

The problem:
Have a look at the images and you'll notice that some detected local minimas (= tree tops) are too close together. There should be only one. This results in multiple watersheds beging generated and I get multiple tree crowns instead of only one.

I have no idea how to treat this. I can use NEAR on the extracted point features to detect those groups, but what next?

One idea would be to detect groups of near-points and only keep the highest (lowest) one and remove the others.

But how?

Or can I filter my initial data more aggressive?


texts: I am currently generating some NDVI imagery using QGIS 2.8.2 and there are a multitude of possible color map scales that can be applied to the imagery. I was wondering if there was a preferred or standard color map to use? USGS, NASA, and others seem to use slightly different scales. I've noticed a wide variety of scales used in academic papers as well. Can any one comment on this?


texts: Can GRASS compute the line graph representation of an input graph?

texts:
I'm interested in computing a line graph for an input graph, preferably in GRASS since I'll be using other GRASS network analysis tools afterwards. I couldn't find any reference to line graphs in the GRASS docs. Is there a different term I should be using in my search or a hidden function that would create such a graph?
For context:

the line graph of an undirected graph G is another graph L(G) that represents the adjacencies between edges of G.

Source: http://en.wikipedia.org/wiki/Line_graph


texts:
I get this question a lot from web developers, and would like to have a short concise answer for them ...

Why do I need ArcGIS Server just to host ArcObjects in a web service?


texts: Exporting classified image from eCognition Essentials Trial license?

texts: How do I export a classified image from the eCognition (Objected based) Essentials Trial?


texts: Error in "Extract values to point" in arcpy

texts:
I am working on spatial interpolation. To extract values for my required points, I am using "Extract values to points" tool from Spatial interpolation toolbox. The tool is working good in ArcMap 10.5. But when I try to use the tool from ArcPy. Same value is generated for each point in test set.
Here is my code:
#point extraction
# Set local variables
inPointFeatures = "test.lyr"
inRaster = "out1"
outPointFeatures = "extractvaluespts.dbf"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute ExtractValuesToPoints
ExtractValuesToPoints(inPointFeatures, inRaster, outPointFeatures,
                  "INTERPOLATE", "VALUE_ONLY")

My code for test point layer is as follow:
   ##XY layer for test
# Set the local variables
in_Table = "test.csv"
x_coords = "LONGITUDE"
y_coords = "LATITUDE"
z_coords = ""
out_Layer = "test"
saved_Layer = r"C:\Users\myc\Desktop\thesis combine\weather parameters\max temp\with elevation\chk1\test.lyr"

# Set the spatial reference
spRef = r"Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"

# Make the XY event layer...
arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer, spRef, z_coords)

# Print the total rows
print(arcpy.GetCount_management(out_Layer))

# Save to a layer file
arcpy.SaveToLayerFile_management(out_Layer, saved_Layer)

I am unable to understand the reason behind this error?
P.S. When I work with geostatistical analyst tool set  and GA to points, it works perfect, but I require to use raster calculator that's why I cannot use Geostatistical analyst. How can I fix this?


texts: Return URL of a query on a feature layer using arcgis API

texts:
I am familiarizing myself with the arcgis API for Python, and I would like to view the query I perform on a feature layer as a URL:

from arcgis.features import FeatureLayer
from arcgis.geometry import filters
# get the data from Boulder, CO and the UDFCD endpoints
city_lims = FeatureLayer("https://maps.bouldercolorado.gov/arcgis/rest/services/plan/CityLimits/MapServer/0")
udfcd = FeatureLayer("https://data.digitaldataservices.com/arcgis/rest/services/UrbanDrainage/UDFCDDataRoutineMaintenance/MapServer/0")

# create the spatial filter test
geom_test = filters.intersects(city_lims.properties.extent)

# perform the query, output to local state plane
cob_udfcd = udfcd.query(out_sr=2876,
                        geometry_filter=geom_test)
cob_udfcd # <-- returns a <FeatureSet>, I want a URL

The documentation mentions additional kwargs that I could enter, and I have a feeling herein lies my answer, but the problem is I have no idea what I would need to enter in order to spit out a URL, and the documentation doesn't offer up much more on kwargs:

  Optional dict. Optional parameters that can be passed to the Query function. This will allow users to pass additional parameters not explicitly implemented on the function. A complete list of functions available is documented on the Query REST API.

How do I tell the ArcGIS API to return a URL which represents a query performed on a feature layer?


texts:
I am new to PostGIS. How can I store 2D elevation data? I have uniform raster and for every [x,y], height is set. I want to be able to use this in my queries (something like SELECT height FROM data WHERE index = [10, 15]). I also need to write this data from my C++ application (it is not earth elevation data, it is some "planet" surface, that will be modified) and data retrieval must be realtime to be used on website.

So far I have come up with 2 ideas:

Use geometry and store 3D data 2)
Create "2D table" with N columns and N rows (not PostGIS, but classic Postgre solution)

Any better ideas or ways how this should be solved?

My data are approx. 3000x1500


texts: Need some hints for a draggable POI-map (Leaflet)

texts:
I try to make such kind of map:
http://www.plotbrowser.com/search.cgi
but first, I can´t see, where they saved the coordinates and the text.
If you look to the sources, there are no coordinates or otehr datas.
I expected that.
But they must be anywhere downloaded, because I can click (already visible/downloaded)
popUps and get new content even when the internet-connecition is blocked
(for sure, offline I can´t drag the map ...).

A good lesson to start is
http://leafletjs.com/reference.html#draggable ??
or other?

I'm working with wampp.

Thanx for tips,
user j


texts: Finding distance in meters between point and path in lat/lon using GeoTools?

texts:
I need to find distance in meters between some point with given latitude and longitude and line (on Earth) defined by several points with lat/lon coordinates.

I read these questions:

Using Linestring distance in meters with GeoTools?
Finding Perpendicular Distance and Minimum Perpendicular Vector between Point and Line using GeoTools and JTS?

In the second it is said that this is impossible with JTS and author decided to use some other utility.

If so, could you recommend any other java library for such problems?

In first, code in the accepted answer doesn't work for me, may be because of that I don't understand which transformation I should choose (and point can be on a big distance from the line, so the transformation may produce not accurate results).

My code for testing code from the first question:

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
    Point spb = geometryFactory.createPoint(new Coordinate(30.3, 59.95));
    Point msk = geometryFactory.createPoint(new Coordinate(37.6167, 55.75));
    Point khabarovsk = geometryFactory.createPoint(new Coordinate(135.0667, 48.4833));

    Point syktyvkar = geometryFactory.createPoint(new Coordinate(50.8167, 61.6667));

    Coordinate[] coords = new Coordinate[] {
            spb.getCoordinate(),
            msk.getCoordinate(),
            khabarovsk.getCoordinate()
    };

    LineString route = geometryFactory.createLineString(coords);
    System.out.println(route.distance(syktyvkar)); //result: 6.881897419667442

    try {
        CoordinateReferenceSystem auto = CRS.decode("AUTO:42001,13.45,52.3");
        MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84,
                auto);
        Geometry g3 = JTS.transform(route, transform);
        Geometry g4 = JTS.transform(syktyvkar, transform);
        System.out.println(g3.distance(g4)); //result: 619842.1091396079
    }
    catch(Exception e) {}
}

The result must be quite close to 170 000 meters.


texts:
I have a raster layer with a corresponding world file. An external program changes the translation parameters in the world file, so I reload the raster by calling

layer.triggerRepaint()

as described in this answer. As desired, the new parameters are loaded, and the image visibly moves to its new position.

However, the canvas clips the raster to the rectangle that it initially occupied, so, for example, if the raster has been shifted to the right by 10 units, all data within 10 units of the right edge are no longer displayed. If I remove the layer and add the raster anew, the whole image is visible in its new location, but this is tedious, and other layer settings (e.g. style, transparency) are lost. I can see in the Python console that the problem is that the layer extent does not change when triggerRepaint is called (whereas the extent of the new layer from re-adding the raster is different from the extent of the original layer). Is there a way to use Python to reset / update the extent of the layer so that the entire raster is displayed in the new location?

I have tried

layer.dataProvider().reloadData()
layer.dataProvider().updateExtents()

to no avail.


texts: SLD style: put all labels at polygon bottom

texts:
I have the timezone polygons in my database and create a layer to it.

Now I try to apply some style but I can't find a way to put all zone labels at bottom.

Already read this http://docs.geoserver.org/stable/en/user/styling/sld-reference/labeling.html and this http://docs.geoserver.org/2.1.1/user/styling/sld-reference/labeling.html.

The result until now:

My style is

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>timezone_polygon</Name>
    <UserStyle>
      <Title>Polygon Timezone Style</Title>
      <Abstract>A Timezone polygon style</Abstract>
      <FeatureTypeStyle>
        <Rule>
          <Name>rule1</Name>
          <Title>Timezone</Title>
          <Abstract>A Timezone polygon style</Abstract>
          <PolygonSymbolizer>
            <Fill>
            <CssParameter name="fill">#f5a6ad</CssParameter>
              <CssParameter name="fill-opacity">0.2</CssParameter>
            </Fill>
            <Stroke>
              <CssParameter name="stroke">#9c0c1f</CssParameter>
              <CssParameter name="stroke-width">0.1</CssParameter>
            </Stroke>
          </PolygonSymbolizer>

       <TextSymbolizer>
         <Label>
           <ogc:PropertyName>zona</ogc:PropertyName>
         </Label>
         <Font>
           <CssParameter name="font-family">Arial</CssParameter>
           <CssParameter name="font-size">9</CssParameter>
           <CssParameter name="font-style">normal</CssParameter>
         </Font>
         <LabelPlacement>
           <PointPlacement>
             <AnchorPoint>
               <AnchorPointX>0.0</AnchorPointX>
               <AnchorPointY>0.0</AnchorPointY>
             </AnchorPoint>
           </PointPlacement>
         </LabelPlacement>
       </TextSymbolizer>

        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: Levelplot(RasterVis) two marginal plot scale?

texts:
How can I add/display scale to two marginal plot?

RAD1998.all <- stack(list.ras)

list.files <- mixedsort(list.files(paste(getwd(), "/1998bil/", sep = ""), full.names = F, pattern = ".asc.bil"))
#list.files

outlist <- substr(list.files, 31, 43)
#outlist

png(paste("1998",outlist[[1]],".png",sep=""))
levelplot(RAD1998.all[[1]])
dev.off()


texts: travel survey multi stop routing problem

texts:
I have a travel diary survey dataset that I am attempting to transform into more useful information.  Currently, the data comes in the format of:

Survey# ,  PlaceVisited#,  TravelMode, TripPurpose, LocationType

1,      1,       NA,        NA,         home

1,      2,       auto,  doctor,     office

1,      3,       auto,  shop,       grocery_store

1,      4,       auto,  NA,     home

where each place the survey respondent visits is a separate record (row) in the survey table, but this is really one multipurpose trip.  What I would like to do with the data is separate each actual trip into a route with various points and purpose information included as (i think) a line using shortest path on a road network (I don't have the route actually taken).  Currently, if a respondent leaves their house, and travels to several different stores by car, each store is entered as a separate record, as is their home starting location.

What I would like to have is one "trip" for each time the respondent leaves their house and comes back, regardless of how many places they visit while preserving the information about each place they visit.  I am wondering what the best way to tackle this problem would be.  I have so far struck out finding prepackaged tools to perform such a task in QGIS, GRASS, ArcGIS, Maptitude, and PostGIS.  Should I approach this as a scripting problem using SQL, scripting a shortest path tool in a GIS, or is there another more efficient way to do it.  Many thanks in advance!


texts:
Trying to calculate aspect for south direction of my DEM. In this manual it is explained how to do it for nothern direction. But when I'm trying to use the expression of this type for southern:

"Aspect@1" >= 135 OR "Aspect@1" <= 225

I'm getting mistake. Both white and black values of resulting raster = 1.

How to change expression to get raster which will give meunderstanding of slopes facing south?
I'm working in QGIS 3.10.4, CRS - Pulkovo 1942 / Gauss-Kruger zone 11.


texts:
I have an hillshade tiff for which the current band type=Bytemay be problematic.

How to change the band type to type=float ?


texts: SpatiaLite: Distance returns NULL when geometries touch

texts:
I use the Python sqlite3 library as well as the SQLite3 CLI.
I am also using the Distance function with the use_ellipsoid parameter, to get the distance between two geometries in meters.

Whenever the geometries touch, I get NULL instead of 0.0 m distance. The exact value of use_ellipsoid does not matter.

An example:

.load 'mod_spatialite.so'
SELECT InitSpatialMetaData(1);
SELECT DISTANCE(MakePoint(13.41, 52.525), MakePoint(13.41, 52.525), 1);

It just returns NULL, respectively, None in Python. If I change one of the points, everything seems to work fine, so it is not a CRS issue.

Is this supposed to be like this? Am I doing or expecting something wrong?


texts:
Due to copyrights I cannot deliver some POIs as plain vector data using GeoJSON in Leaflet or OpenLayers. I'm therefor looking for some kind of sane obfuscation, to make it harder for normal users to get to the coordinates.

Creating a new map projection might not be a good solution, since the webmapping library has to contain the code to reproject to the map's CRS.
MVT maybe a better solution. Of course the original coordinates still could be reversed, albeit with greater effort.
Seems like the most secure solution is to precalculate the POIs as raster tiles with a WMS(T) (of course with the map icons). Interactivity would come via UTFgrid.
I also was thinking about SVG or local leaflet x-y coordinates. Here reverse engineering is still possible, but the effort increases extremely.

Anyone with a better idea?


texts: Data Transfer not working, Files to Pathfinder another way?

texts: Typically I would use the Data Transfer Utility to retrieve files from my Trimble Geo7x and run them through GPS Pathfinder to create usable shapefiles. Currently I cannot get the transfer utility to talk to Terra Sync resulting in my inability to load files into GPS Pathfinder. I am working on fixing that but in the meantime, because I am on a tight schedule, does anyone know if there is an alternate way to retrieve files from the device and put them directly into GPS Pathfinder without the transfer utility? I need to be able to correct and create shapefiles for ArcGIS as soon as possible.


texts:
I'm trying to convert a 20 GB LZW compressed geotiff to a tilepyramid using gdal2tiles.
Unfortunately, gdal2tiles seems not to be able to handle this, after running it for 10 hours it still hasn't reached the first "dot" yet.  I tried this with a smaller (1GB) file and it worked fine.  My question now is, is there any other software that allows me to do this (commercial or open source)?  One can create tiles from GlobalMapper, but I cannot load the 20GB file into memory.


texts:
GWR:
What is the difference between "Local R2" and R2 is the result of GWR. In the output table, I have an R2=%70 but in another table, the maximum Local R2 has been reported as %21.


texts: Count the polygons of the raster given by the kmeans algorithm

texts:
I am struggling with a Python problem. My plan is to compute the clustering of the kmeans algorithm for an image and to create polygons from this classification. My goal is to get the total number of Polygons and the area they are capturing.
The computation of the kmeans classification is working without any problem, but I have problems to create Polygons for a selected cluster. My code looks like this:
###import the necessary libraries
import cv2
from sklearn.cluster import KMeans
import numpy as np
from rasterio import features
from shapely.geometry import shape

###import the image, you can test it with any image
img = cv2.imread("Test.jpeg")
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

###image preperation
x, y, z = img.shape
im = img.reshape(x*y, z)
im = np.float32(im)

#final cluster computation
kmeans = KMeans(n_clusters=10, n_init=2)
kmeans = kmeans.fit(im)
means = np.float64(kmeans.labels_)

#select one cluster
masked_data = np.ma.masked_where(means != 0, means)

# try to compute polygons based on the classification
m=masked_data.data.reshape((x,y,1))
m=m.astype(np.float32)
mypoly = []

for vec in rasterio.features.shapes(m):
    mypoly.append(shape(vec[0]))

When I export the mypoly array, I get a strange looking polygon. I exported it like this:
import geopandas as gpd
my_gdf = gpd.GeoDataFrame( geometry=mypoly)
my_gdf.to_file("Example.shp", driver='ESRI Shapefile')

Any idea how to fix this?


texts: Determining percent of polygons within other polygons using ArcGIS Pro?

texts:
I am using ArcGIS Pro.

I have two layers 'Habitat Quality' and 'Connectivity'. I am trying to find out the percent of habitat quality polygons within the connectivity polygons. The connectivity layer is already a joined layer and so its attribute table has a lot of nulls because the layers were not perfectly matched.


texts:
I like the idea of Virtual Layers, but when they are too big and constructed from complex queries, QGIS takes a lot of time to redraw it (it seems like it recalculates the entire layer at every zoom/pan/etc).

My workaround is to export Virtual Layer to GeoPackage format and continue to work with that.

In my case, Virtual Layer is often just an intermediate to be processed with further tools and I do not need it long term. I would rather save it as "Temporary scratch layer", just like many Processing tools do.

Is there a way to do so?


texts: Joining a row of data to multiple polygons (zip codes)

texts:
I received some data from a partner and was trying to do a visualization. I'm using ArcGIS / ArcMap 10.2

I have a table where each row (first column is the ID) represents a facility and all the US zip codes that intersect within a certain radius of that facility.

The zip codes are listed in a concatenated field, called "RSAComponents"

What could I do to join shapefiles representing multiple associated zip codes to each row?

I suspect I could do some iteration via ModelBuilder, but I'm still a novice at using it.


texts: Buffer not providing accurate geodesic buffers in ArcGIS Desktop 10.0?

texts:
This appears to only be a problem when using ArcGIS Desktop 10.0.

I created a shapefile following these criteria as recommended for Performing Geodesic Buffering:

The input feature class must be a point or multipoint dataset. (shapefile is a point)
The input feature class must have a geographic coordinate system. (GCS: Longitude/Latitude (WGS 84))
The buffer distance must be specified as a linear unit such as kilometers, miles, or feet. (Radius field (String) added with inputs 3 kilometers, 4.5 kilometers...)

I am not, however, getting accurate geodesic buffers. If I perform a buffer with the same shapefile and input a Linear Unit (3 Kilometers) instead of specifying by the field, I get accurate geodesic buffers.

This was the buffer used with a specified distance which gave me accurate geodesic buffers:

arcpy.Buffer_analysis("test","C:/Temp/testBuffer.shp","3 Kilometers","FULL","ROUND","LIST","Site_Name")

This was the buffer used with a specified field and is NOT giving me accurate geodesic buffers:

arcpy.Buffer_analysis("test","C:/Temp/testBuffer2.shp","Radius","FULL","ROUND","LIST","Site_Name")


texts: Clustering point data to obtain maximum count groups using QGIS?

texts:
I am working with a data set of over 3000 points which I need to cluster/group based on a distance criteria.

What I would like to do is build clusters from the point data which contain the maximum number of points possible to group in a 50 m radius. To be clear I don't want to group any points that are within 50 m of each other.

The way I envisage it is taking a 50 m radius polygon and centring it one by one over each point and counting how many other points fall inside the polygon. At a simple level this would return a value for each point saying how many points are within 50 m of it. I could then choose groups based on which points return the most other points with 50 m.

I would however like to be even more iterative than this as doing the above does exclude a point falling in multiple 50 m polygons. I would like a process which looks over the data and provides the most inclusive 50 m clusters for the data set as a whole.

I'm not sure if QGIS has a tool which could do this at all or whether I should look to develop a process in Excel or elsewhere.


texts: Calculating pixels coordinates from lat/lon

texts:
I'm trying to develop a WPS that would print maps to images and PDFs. I don't want to use Mapfish because of its limitations, and because of some customization I need to add to the printed maps.

My method goes something like the following:

1- User selects a BBox to print.

2- The WPS converts the BBox coordinates to pixels coordinates.

3- Then calculates the resulting image's width and height using the calculated pixels coordinates.

4- Then call a WMS GetMap request to get the map image.

The problem is that the resulting image is stretched longitudinally.

I have reviewed the code many times and it seems that the problem is in the calculations of the pixels coordinates, it returns incorrect values which result in incorrect width and height, which in turn makes the image stretched.

Here is the code I use to convert from latitude/longitude to pixels coordinates, I took it from an opensource project called MOBAC.

Note  that my layer is in EPSG:900913.

private void init() {
    for (int zoom = 0; zoom < worldSize.length; zoom++)
        worldSize[zoom] = tileSize * (1 << zoom);
}

/**
 * Returns the absolute number of pixels in y or x, defined as: 2zoom * tileSize
 *
 * @param zoom
 *            [0..22] (for tileSize = 256)
 * @return
 */
private int getMaxPixels(int zoom) {
    return worldSize[zoom];
}

private int falseNorthing(int aZoomlevel) {
    return (-1 * getMaxPixels(aZoomlevel) / 2);
}

private int cLatToY(double lat, int zoom) {
    lat = Math.max(MIN_LAT, Math.min(MAX_LAT, lat));
    lat = Math.toRadians(lat);
    lat = -Math.log(tsfn(lat, Math.sin(lat)));
    int mp = getMaxPixels(zoom);
    int y = (-1) * (int) (mp * lat / (2 * Math.PI));
    y = y - falseNorthing(zoom) - (y > 0 ? -1 : 1);
    y = Math.min(y, mp - 1);
    return y;
}

private int cLonToX(double lon, int zoom) {
    int mp = getMaxPixels(zoom);
    int x = (int) ((mp * (lon + 180l)) / 360l);
    x = Math.min(x, mp - 1);
    return x;
}


texts:
This is the example code to load a Landsat 7 image
img = ee.Image('LANDSAT/LE07/C01/T1_SR/LE07_137044_20200813')

print(img.getInfo())

I got this error:

HttpError                                 Traceback (most recent call last) ~\Anaconda3\envs\gee\lib\site-packages\ee\data.py in
_execute_cloud_call(call, num_retries)
    344   try:
--> 345     return call.execute(num_retries=num_retries)
    346   except googleapiclient.errors.HttpError as e:

~\Anaconda3\envs\gee\lib\site-packages\googleapiclient\_helpers.py in positional_wrapper(*args, **kwargs)
    133                     logger.warning(message)
--> 134             return wrapped(*args, **kwargs)
    135

~\Anaconda3\envs\gee\lib\site-packages\googleapiclient\http.py in execute(self, http, num_retries)
    906         if resp.status >= 300:
--> 907             raise HttpError(resp, content, uri=self.uri)
    908         return self.postproc(resp, content)

HttpError: <HttpError 400 when requesting https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/value:compute?prettyPrint=false&alt=json returned "Image.load: Image asset 'LANDSAT/LE07/C01/T1_SR/LE07_137044_20200813' not found.">

During handling of the above exception, another exception occurred:

EEException                               Traceback (most recent call last) <ipython-input-31-8676335bc0d4> in <module>
      1 # Print image object WITH call to getInfo(); prints image metadata.
----> 2 print(img.getInfo())

~\Anaconda3\envs\gee\lib\site-packages\ee\image.py in getInfo(self)
    116           properties - Dictionary containing the image's metadata properties.
    117     """
--> 118     return super(Image, self).getInfo()
    119
    120   def getMapId(self, vis_params=None):

~\Anaconda3\envs\gee\lib\site-packages\ee\computedobject.py in getInfo(self)
     93       The object can evaluate to anything.
     94     """
---> 95     return data.computeValue(self)
     96
     97   def encode(self, encoder):

~\Anaconda3\envs\gee\lib\site-packages\ee\data.py in computeValue(obj)
    708             body={'expression': serializer.encode(obj, for_cloud_api=True)},
    709             project=_get_projects_path(),
--> 710             prettyPrint=False))['result']
    711   return send_('/value', ({'json': obj.serialize(), 'json_format': 'v2'}))
    712

~\Anaconda3\envs\gee\lib\site-packages\ee\data.py in
_execute_cloud_call(call, num_retries)
    345     return call.execute(num_retries=num_retries)
    346   except googleapiclient.errors.HttpError as e:
--> 347     raise _translate_cloud_exception(e)
    348
    349

EEException: Image.load: Image asset 'LANDSAT/LE07/C01/T1_SR/LE07_137044_20200813' not found.

However, when I tried with example data from https://developers.google.com/earth-engine/python_install#printing-objects it works and prints the metadata information. What might be the problem here?


texts: Choosing correct zoom level to display given area in square km using Leaflet?

texts:
I would like to set the zoom on a Leaflet map based on an area (in km^2).

I am looking for a simple function that would approximately convert an area of
510,064,472 km to zoom level 0 and (if I read the table of zoom level at https://wiki.openstreetmap.org/wiki/Zoom_levels correctly) an area of 12,321 to zoom level of about 8.5 and so on.


texts: Creating Contours from a Classified Raster - ArcGIS

texts:
I have a group of classified raster files, showing point densities at 3 intervals. I want to create contour lines based on this data at these same intervals. I'm guessing this should be fairly easy, but how? I'm on ArcGIS 10.1.

EDIT:
I have about 100 raster files, each of them is a group of points which I've performed a Kernel Density Estimation on, so is now a raster. I've classified these into 3 classes 15-50, 50-100, and 100-250 (there is also a 0-15, but I don't need this included). How would I go about creating contours at each of the limits of the 3 classes I want to keep?
Martin - Deleting extra contours would be fairly time-consuming. The operation would preferably be done as a batch - they all have the same class symbology and they would all have the same contour limits too.


texts: Accessing MySQL data from python built-in console of QGIS?

texts:
I'm writing a QGIS Plugin, and somewhere I want to present data from a MySQL database. It's simple data(strings,etc), not Geospatial, or vector data. I'm familiar with querying databases from python scripts, and now I want to present data inside a TextBrowser QT Widget inside QGIS.

However, I can't import MySQLdb from built-in console, neither from my xxx.py file inside my plugin folder. From the Apple Terminal i can import MySQLdb with no problems.

I've read Why can't I import my script from the python console?, but still have no success, because when I call import MySQLdb inside a script that resides on a folder added to sys.path, it still fails.

Calling  import MySQLdb from built-in console tells:

import MySQLdb
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "build/bdist.macosx-10.8-intel/egg/MySQLdb/__init__.py", line 19, in <module>
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "build/bdist.macosx-10.8-intel/egg/_mysql.py", line 7, in <module>
  File "build/bdist.macosx-10.8-intel/egg/_mysql.py", line 6, in __bootstrap__
ImportError: dlopen(/Users/Arthur/.python-eggs/MySQL_python-1.2.3-py2.7-macosx-10.8-intel.egg-tmp/_mysql.so, 2): no suitable image found.  Did find:
    /Users/Arthur/.python-eggs/MySQL_python-1.2.3-py2.7-macosx-10.8-intel.egg-tmp/_mysql.so: mach-o, but wrong architecture

ACKNOWLEDGE:

..mach-o, but wrong architecture -> this reveals a mismatch between python-mysqldb-mysql architecture!

The problem is that MySQLdb, does not support python2.7, as compiling for x86_64!
Even if python and mysql-server are 64bits. Switching to python2.6, MySQLdb compiles a x86_64 binary; but QGIS 1.8 is built against python 2.7, so still does not work.
Is it possible to select python framework inside QGIS?

This was very helpful -> https://stackoverflow.com/questions/12327737/error-importing-mysqldb-in-python-2-7-2-on-mac-osx-10-8-1


texts: Split line at point doesn't split lines less than 30 meters

texts:
I have a bunch of lines in line layer that should be split at midpoint.
I tried to use both Feature Vertices to Points and Feature to Point to get a new layer with point at the middle of lines. I even used topology check to fit points and lines at coordinate grid.

Eventually lines that have a length more than 25-30 meters are split while others are not! And I need to split all the lines in layer.

On a picture above you can see that two last lines at right are not split.

I also tried to use different projections but it didn't help.
Is there any solution or this is a bug of a tool whcih cannot split lines with length that is less than 30 meters?

UPDATE

Found out that export lines to shapefile gives positive result, i.e. lines can be split in that way.
But I try to use only GDB file format so that and the final result of this operation will be presented in ModelBuilder.


texts: Extracting Multi Values to Points temporally as well as spatially from NetCDF raster layer using ArcGIS Desktop?

texts:
I have created a map in ArcGIS 10.6 which includes point data from ARGOS satellite tagged turtles. Each point has a date/time field in a dd/mm/yyyy hh:mm format. I also have a sea surface temperature (SST) as a NetCDF from ERA-Interim with 6-hourly resolution which I have added to the map using 'Make NetCDF Raster Layer'.

Is there a way in ArcGIS to 'Extract Multi Values to Points' temporally as well as spatially from this data?

The time in both datasets will not match exactly so, I need to extract the SST data which is closest to the point time.


texts: Start a loop from information in a feature property

texts:
I have a feature collection of polygons that contains a property ("yrbegun") indicating the year I want to use to do some manipulation. I would like to start a loop from the year before to the year after, i.e. if the property indicates 2005, I want to loop from 2004 to 2006

// Determine the year the project was started
var projectyear = ee.Number(inFeature.get("yrbegun"))

// Define the "boundaries" of the loop
var startyear = projectyear.subtract(1)
var endyear = projectyear.add(1)

// Run the loop
for (var yr= projectyear; yr <= endyear; yr = yr + 1) {

  print('This is year ' + yr)
  print(yr)

}

However, the loop never runs. I tried to define projectyear without the ee.Number() function, but it also does not work.


texts: Split polygons by polyline but with specified attribute

texts:
I have a feature class of polygons (ImportPoly), each of them is intersected by line from (ImportLine) and then new polygons created in (OutPoly) feature class with all attributes saved. I need to cut each polygon with only one line that has same [cut_field] as [cut_field] of polygon. Methods described here and here don't really work for me, because each poly must be splitted only with corresponding line.
Is there a way to write something like "search cursor in ImportPoly and ImportLine, if cut_fields are equal then use Feature to Polygon tool for them"?

Edit: solution by @Hornbydd works perfectly in any way but time. Each iteration took roughly 10 seconds, which would become great amount of time for large datasets.


texts: MapGuide Maestro - Input string was not in a correct format

texts:
I facing an issue by using MapGuide Maestro.
I installed MapGuideOpenSource 2.6.0.8316 x64 and MapGuideMaestro Release 6.0a3.
When I started MapGuide Maestro a first windows info appear and say:

I create my folder and I add my shapefile resource, double click on my layer and in order to modifed my feature representation and labels.

When I click on previous resource icon, an error message show up.

Somebody have an idea?

I found several tutorial online and look around for a solution but nothing relevant.


texts: Google Earth Engine: Working with imagery with "no bands" --how to select or rename a band

texts:
I am trying to perform various operations on MODIS datasets in google earth engine.
Every operation I try to perform is prevented by the requirement in code to select a band to work with.
This collection: MODIS/MCD43A4_006_NDVI'
has no bands. It shows up in the console as
bands: [] when I perform any operation that requires selecting or manipulating bands.
part of the code I've tried is this:
  .filter(ee.Filter.calendarRange(startyear,endyear,'year'))
  .filter(ee.Filter.calendarRange(startmonth,endmonth,'month'))
  .sort('system:time_start')
  // .filterBounds(ROI)
  .select('NDVI');

I can't select NDVI, which I need later in the code, because I can't select anything with this dataset--there are no bands. I can't add a band, since the information is already there, NDVI and day of year, but apparently not stored in a band.
How can I use this dataset in operations that require selection of a band?


texts: Killing the correct ArcMap process from Windows task manager

texts:
I have 3 ArcMaps opened, one of them has a process which is taking a long time but still running. Another one hangs and I can't close it.

I want to kill the hung process, but how do I know from the task manager which one it is? I don't want to kill the one which has the long running process.


texts: How to remove first two and last two characters from each label in ArcMap

texts:
My too long labels take up too much space in regard to the small counties polygons.
I want to remove first two and the last two characters from each single label.
Currently the labels on the map look like this:
PL1607000 PL1816000 LV801200 LT029900 EE700000 etc.

Because in each single label the first two characters are two-letter country code, and the last two characters always two zeros, to make map more clean I want my labels like this:
16070 18160 8012 0299 7000 etc.

I use Maplex Label Engine. My label field is CTYCode, What expression should I type in Label expression dialog box?
Is it possible to also make stacked labels from this labels? I want to the first two characters were at the top, and the rest at the bottom. Like this:

16       18     80    02   70
070    160    12    99   00

  etc.

How do I write two different expressions. With and without stacking labels, I do not know about programming at all and I would not know how to change the expression.


texts: Query a region within a larger region with GeoDjango REST API

texts:
I have a working Django REST API which can query all the counties in a country.
//My model
class kenya_counties(models.Model):
  objectid = models.BigIntegerField()
  id_field = models.BigIntegerField()
  county_nam = models.CharField(max_length=80)
  const_code = models.BigIntegerField()
  constituen = models.CharField(max_length=80)
  county_cod = models.BigIntegerField()
  shape_leng = models.FloatField()
  shape_area = models.FloatField()
  geom = models.MultiPolygonField(srid=4326)

  def __str__(self):
      return self.county_nam

//My Serializer
class CountiesSerializer(serializers.ModelSerializer):
  class Meta:
     model = kenya_counties
     fields = ('objectid', 'id_field', 'county_nam', 'const_code', 'constituen', 'county_nam', 'county_cod', 'shape_leng', 'shape_area')

//My view
@api_view(['GET', 'POST'])
def counties(request):
   if request.method == 'GET':
       data = kenya_counties.objects.all()

       serializer = CountiesSerializer(data, context={'request': request}, many=True)

       return Response(serializer.data)

   elif request.method == 'POST':
       serializer = CountiesSerializer(data=request.data)
       if serializer.is_valid():
           serializer.save()
           return Response(status=status.HTTP_201_CREATED)

       return Response(serializer.errors, status=status.HTTP_400_BAD_REQUEST)

//My url
    re_path(r'^api/gis/counties/$', views.counties),

I get .
How can I query just one county and display it on a react-leaflet map as shown below?
<MapContainer center={position} zoom={13} scrollWheelZoom={false}>
    <TileLayer
      attribution='&copy; OpenStreetMap contributors'
      url="https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png"
    />
</MapContainer>


texts: Radius of Curvature Python/Arc, how access x,y coordinates

texts:
I have a polyline in ArcMap 10 and I am able to script the length of the line in Python using this bit of code:

#Calculate Radius of Curvature
#try:
# Import system modules
import sys, string, os, arcpy, math, shutil

#Set Environment Workspace
arcpy.env.workspace = "C:\data\ColRiver\ColRiver.gdb"

# Create search cursor
rows = arcpy.SearchCursor("c:\data\ColRiver\Colriver.gdb\ColumbiaRiverPolygon")

# Calculate the total length along the flowline
length = 0

shapeName = arcpy.Describe("ColumbiaRiverPolygon").shapeFieldName

for row in rows:
    # Create the geometry object
    feat = row.getValue(shapeName)
    length = length + feat.length

# Print the total length of flowline
print length

Now, I know that to get the radius of curvature for the line I will have to divide the total length by the straight line distance of the start point and end point.  I understand that the way to get distance would be:

((x2-x1)^2 + (y2-y1)^2)^0.5

What I am unsure of is how to automate getting the start point and end point coordinates of x1, y1 and x2, y2 in order to get this straight line or euclidean distance. I plan to implement part of the code in a tool/script for ArcMap and don't want to have to type in the coordinates every single time. Any suggestions for a Python newbie?


texts: Is there a way to copy the mxd file with all of the files in that in ArcGIS?

texts: Obviously the mxd doesn't contain the gdb files, and the linked raster files. So when I copy a project, I need to search manually the 49 raster maps from the other 4000. Is there an easier way to copy this?


texts: Creating a domain based off specific attribute selection or lookup table in ArcGIS/Access?

texts:
Is it possible to have a domain on a database in ArcGIS that creates a drop down pick list on an attribute then have the next field have a domain based on that section?

For example:

Field A - choose between cat, dog, cow, chicken, snake
Field B - choose between white, black, red, brown ONLY IF I PICK 'CAT'
Field B (may also) - choose between purple, yellow, blue ONLY IF I PICK 'DOG'
FIELD C - stripes, spots (ONLY IF I PICK 'CAT' THEN 'WHITE'

and so on........

I do not want to create a script and run that   to update the attribute table in arc every time a selection of list of attributes is made.


texts: Joining more than two lines that intersected in one point using FME

texts:
i have point connects with four lines  in dgn file
Each line has his unique id ,iam usng FME 2018 software

i used line combiner for the map ,i found only ,each two lines is able to connect together .but when more than two line are intersecting in one point , line combiner will not able to connect them .

i tried a lot for some days but it does not work .

Do you have a suggestion for how I can do that using FME?


texts: OSMasterMap.ttf not recognised in QGIS - macOS Sierra

texts:
I am currently learning how to work with OSMastermap Data with postgreSQL and postGIS. After much tinkering I have finally got OS Translator II working and have the Topographic data converted and in the database, with style applied.

The only issue I have now is I can't seem to get OSMasterMap.ttf to work properly within QGIS. It recognises that the font is present, but it does not apply the font when rendering, it reverts back to normal letters.

When I go to select a font it can be selected from the drop down, but the symbols are not used.

The font is installed on the system, as it can be seen in the font book.

Does anyone know if there is a solution to this problem?


texts: QGIS panel will not dock (or too small and will not undock)?

texts:
I detached my Layers panel from QGIS (version 3.6.2) and now it will not dock.

I have already tried double clicking the top bar, this works for every window/panel apart from Layers. It simply will not go back.

Additionally if you look below the toolbar, this is the Vector Toolbar. I cannot get it to undock or resize to even show me the icons, they are lost forever.

Is there a config file or setting I can delete to revert to defaults for this?


texts: Extracting Raster value from sampling point

texts:
I have Sample RGB combination from 3 different Bands

// Input RGB imagery
var image = ee.Image("users/k2yulianto/Exercise/Ras/Imagery");

// Use Tree bands from RGB Imagery.
var bands = ['b1', 'b2', 'b3'];
var b1 = ['b1'];
var b2 = ['b2'];
var b3 = ['b3'];

// Display the inputs imagery
Map.addLayer (image, {bands: ['b1', 'b2', 'b3'], max: 1, gamma:1, opacity:0.80} , 'RGB imagery',true);
Map.setCenter (101.4476, 0.4925, 12);

and have several sampling point on top of the imagery :
// Sampling point on top of RGB imagery
var Sampling = ee.FeatureCollection("users/k2yulianto/Exercise/Sampling/Sampling");
Map.addLayer (Sampling, {opacity:0.80},'Sampling point', true);

i want to extract raster value in the each bands from imagery using my sampling point
and exported in to shapefile /csv.
the shapefile/csv should be contain column of original bands combination and each band value in it:

ID | Latitude | Longitude | b1 | b2 | b3

// Empty Collection to fill
var ft = ee.FeatureCollection(ee.List([]));

var fill = function(img, ini) {

// gets the values for the sampling points in the current image
var getvalue = image.reduceRegions(Sampling, ee.Reducer.first(),30);

};

// Iterates over the image
var newft = image.iterate(fill, ft);

// Export
Export.table.toDrive(newft,
"anyDescription",
"anyFolder",
"anyNameYouWant");

When i run my scripts and i got error warning

"image.iterate is not a function".

How can i solved this problem.
i need to export the value raster to shapefile or csv
that the column will contain :

ID | Latitude | Longitude | b1 | b2 | b3

here is the link of my original scripts:

https://code.earthengine.google.com/6e4200dfe8c8b130430f8e4367712ebf


texts:
I was taught that you need to 'Analyze' each dataset in enterprise databases as much as you can. It was one of those things that I just did, and which I never thought much about. However, I recently looked up the tool I always use, to see what it actually does: ArcGIS Pro Analyze.
It states:

Updates database statistics of business tables, feature tables, and
delta tables, along with the statistics of those tables' indexes.

Does this really impact any performance in a measurable way?
Is this something that really is needed?
It seems to me that it is something that could be done (and slightly improve database performance), but really does not do a lot.


texts: Earth Engine: different results with "covariance" and "centeredCovariance" on de-meaned data?

texts:
I am trying to compute a covariance matrix in earth engine. EE offers two reducers:

ee.Reducer.centeredCovariance
ee.Reducer.covariance.

My understanding is that ee.Reducer.covariance on the raw data should be the same as ee.Reducer.centeredCovariance on a band-demeaned data (from help file: this reducer requires that the data has been mean centered.) This is also suggested in this EE script, which is mentioned on the official documentation (eigen anaysis). However, computing both approaches in the same script leads to quite different results!?
How to explain these different results? Is my expectation that the two should the same wrong, did I make a mistake in the code, or is there a problem in the way ee.Reducer.covariance is computed (using apparently this algorithm)? Note that there might be a difference in the degrees of freedom adjustment (1/N vs 1/(N-1)) but with 1896290 pixels in the zone, this difference ratio is at most 1.000001, while difference ratio in the results is easily 1.122402!
Code (link):
// Script mainly copy-pasted from https://code.earthengine.google.com/f80ead7606e2feeaa3744cf3d03d761d,
//which is mentioned on https://developers.google.com/earth-engine/guides/arrays_eigen_analysis

// Load a landsat 8 image, select the bands of interest.
var image = ee.Image('LANDSAT/LC8_L1T/LC80440342014077LGN00')
  .select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11']);

// Display the input imagery and the region in which to do the PCA.
var region = ee.Geometry.Polygon([[[-122.77, 37.70], [-122.77, 37.44],[-122.10, 37.44],[-122.10, 37.70]]], null, false);
Map.centerObject(region, 10);
Map.addLayer(ee.Image().paint(region, 0, 2), {}, 'Region');
Map.addLayer(image.clip(region), {bands: ['B5', 'B4', 'B2'], min: 0, max: 20000}, 'Original Image');

// Set some information about the input to be used later.
var scale = 30;
var bandNames = image.bandNames();

// Mean center the data
var meanDict = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
});

var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);

// Compute the covariance of the bands within the region.
var covar_cent = centered.toArray().reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
  });

var covar = image.toArray().reduceRegion({
    reducer: ee.Reducer.covariance(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
  });

print(ee.Array(covar_cent.get('array')).slice(1, 0, 1))
print(ee.Array(covar.get('array')).slice(1, 0, 1))


texts: Checking intersection of two shapefiles using ArcPy

texts:
I am new to Python programming.
How do I check the intersection of two shapefiles using ArcPy?
If the shapefiles intersect, it should return True else False.
Below is the code that I have written :
import arcpy

arcpy.env.workspace = "D:/Qasim/mn_corcoran"

arcpy.Intersect_analysis(r'D:/Qasim/mn_hopkins/abc.shp #;D:/Qasim/mn_corcoran/xyz.shp #', r'D:/Qasim/mn_hopkins/asd.shp', 'ALL', '#', 'INPUT')

if arcpy.Intersect_analysis == True:

 print("True")

else: arcpy.Intersect_analysis == False

print("False")


texts:
I am trying to load a DEM (img file + prj file) into QGIS 3.12.1 and use it as a layer over OSM. However I can't seem to get the CRS right, the OTF-conversion does not show the layer on the map.

I am using EPSG:3857 as project default and the read CRS from the prj file with OTF CRS conversion enabled:

PROJCS[
    "Transverse_Mercator",
    GEOGCS["GCS_Geographic Coordinate System",DATUM["D_MGI",SPHEROID["Bessel_1841",6377397.155,299.1528131060786]],
    PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],PARAMETER["scale_factor",1],
    PARAMETER["central_meridian",10.33333333333333],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-5000000],
    UNIT["Meter",1]
]

Does anybody have an idea how I can transform this CRS to EPSG:3857?

Here are the files I am trying to load: https://vogis.cnv.at/geodaten/?service=files&t=418171be46e40175bfe171edfc4097dd&path=%2FGelaendemodelle%2FHoehenmodelle%2FVlbg_Umgebung

Thank you


texts: Configure the Search Function in QGIS Web Client

texts:
I have set up the QGIS Web Client but do not know how to get the search function to work. I read the instructions in the guide but can't really understand what it is saying:

 6.1.2. PHP configuration file
    =============================

Configuration for the services is stored in  `config.php`.

Example:

  /****************************
   * Map rewrite configuration
   */
  // Prefix map name with path
  #define('MAP_PATH_REWRITE', '/home/xxx/public_html/QGIS-Web-Client/projects/');
  // Append .qgs to the map name
  #define('MAP_PATH_APPEND_QGS', true);

  /**************************************
   * search configuration
   */
  // Configuration for searchable layers
  $searchlayers_config = array(
      // Key is layer name
      'Country' => array(
          // SQL for text search: where to search
          'search_column' => 'name'
      )
  );

  // Default search tables
  define('DEFAULT_SEARCH_LAYERS', 'Country');
  // Limit search results
  define('SEARCH_LIMIT', 100);

  /**************************************
   *  Get legend configuration
   */
  // Cache expiry time in seconds 0=never cache
  define('GET_LEGEND_CACHE_EXPIRY', 60*60);
  // Cache directory, defaults to dirname(__FILE__) . '/legend_cache'
  define('GET_LEGEND_CACHE_DIRECTORY', null);
  // Defaults to current URL + '../cgi-bin/qgis_mapserv.fcgi?'
  define('WMS_ONLINE_RESOURCE', null);

  /* End configuration */

QGIS Web Client needs to know where to find the scripts, since most
configuration is read from the project file, this must be passed in the
query string, the file where this parameters are set is
GlobalOptions.js see the example below:

  // Adds project_map, read value from query string
  var project_map = Ext.urlDecode(window.location.search.substring(1)).map;

  var searchBoxQueryURL = '../php/search.php?map=' + project_map;
  var searchBoxGetGeomURL = '../php/search_geom.php?map=' + project_map;

Can anyone enlighten me what must I do, which part of the code to edit? Thanks.


texts: How to remove lines made with a PostGIS query automatically if origin is not in the map window?

texts:
I have a PostGIS point layer with labels. I moved these points because they were too close.

To connect label and point, I created a line vector layer which was made with a query in the DB-Manager (see below).

If a point is not visible because it's outside of the QGIS window, but the label is still inside, the label is not visible anymore. But, unfortunately, these lines are still visible. What do I have to add to the query to enable these lines if their origin is not visible?

Query:

SELECT
  id,
  ST_Makeline(St_setSRID((geom),2056), St_setSRID(St_Point("X"::numeric, "Y"::numeric),2056)) as "geom"
FROM
"database"."table"
WHERE
  "X" IS NOT NULL AND
  "Y" IS NOT NULL


texts: I have placed 10 geopackage files locally on GeoServer. I also added a plugin for vector tiles. Tell me the code to open them via Leaflet and that when you click on the object there is information.


texts: Adding segments3d to a view 3D from lidr

texts:
I want to add segments represent the tree of apex tree, how to add segments3d to the plot below ?

library(lidR)
library(rgl)

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyz", filter = "-drop_z_below 0")

ttops <- tree_detection(las, lmf(ws = 5))

plot(las) %>%
  add_treetops3d(ttops) %>%
  segments3d(...)

What are arguments of ... ?


texts:
I have been fighting with this for a week or so.  I  believe I have viewed the related questions on stack exchange, tried some of the suggested solutions that seemed appropriate and have found none that seem to address the same issue I seem to be having. So, I do not think this is truly a duplicate question.

Caveats: New to QGIS, Not a programmer, very much in the learning stages with all elements of QGIS.  Using the program for fictional world creation.  Have had considerable success but realize I am trying to use the program for a purpose that is not its main function.

Workflow to Date

I have downloaded a section of topographic data from Open Topography.
I converted the data to an excel format.
I have altered the x and y coordinates by a constant (basically taking a section of Peru and relocating it to 140 degrees East longitude and 40 degrees North latitude in a 5 degree section (i.e. 140-135E, 40-35N)
I saved the file as csv.  I used the delimited text file layer function and successfully imported it into QGS 3.10 - both with and without GRASS 7.6.  It appears (in the correct latitude and longitde I wanted, but with the original data for elevation) as a layer of very dense points (see photos.)
I have extracted vertices, extracted coordinates, made a successful contour map and seen the attribute table fully populated with almost 900,000 data points with x y and z coordinates where they should be.

The Problem

I cannot create a useful DEM raster from the data.  In the past I have used drawn polygons and done the same extractions a above and then done TIN interpolations, with success, to create DEMs.  TIN does not work in this case, and the row x column generated is 3 x 4, when it should be in the tens of thousands.

Attempted work arounds

I have saved the layer as a shapefile and as a geopackage.
I have tried extracting the vertices and then added the geometry attributes,
I have tried adding an elevation field and setting the z coordinate to that (this solved an earlier problem.)
I have converted the file to contours and extracted points from that again, in an effort to reduce the data.
I have attempted TIN and IDW interpolations - but with the density of the points, neither seems appropriate somehow? - too close together for TIN and too dense to need IDW?
I have tried various GRASS and other transitions (v.surf.idw, r.surf.idw, vector to raster in the rasterize menu) but none of these have worked in terms of changing that low numbers of rows and columns, which I (possibly erroneously) assume means that it is somehow not reading the data?

Screenshot

I have included a screen shot of the data points as they appear when I load the file.  I have done them at a large (1:48,000 and small 1:1137) scale.

The Question(s):

What data inputs would you need to see to help me solve this issue?
Is there anything you might recommend I do differently here?
While I have been looking for days at tutorials that don't seem to give me an answer, are there any tutorials you would recommend?

Additional Photos

I am adding several screenshots

The CSV load screen - you can see the organization of the table here

The attribute table for the layer once it is loaded

The settings as I save it as a new layer

The TIN interpolation settings

The raster disaster :)

A final shot showing the contour layer that can be generated from the CSV data.

Raster disaster number 2 after I tried setting everything to UTM and them the columns to 400 manually...

Reloaded the file, setting everything to UTM zone 7N - top bit looks good - will try it again


texts: zoom to cluster onClick & have hover event

texts:
I've made a simple mock-up trying to combine two openlayers examples. I want to have clusters on which I can hover (so I can display an infowindow about it's content) and I want users to be able to click on them to zoom to these clusters.

The last part I've allready made. Take a look at this jsFiddle
http://jsfiddle.net/BzcAD/2/

You might recognize the following two examples:
http://www.openlayers.org/dev/examples/highlight-feature.html
 http://www.openlayers.org/dev/examples/strategy-cluster-threshold.html

Now the jsFiddle halts on the error : "TypeError: this.layer.getZIndex is not a function"

You can stop this error by uncommenting the lines regarding the selectCtrl. Like so:

// var selectCtrl = new OpenLayers.Control.SelectFeature(
//     {clickout: true}
// );
   map.addControl(highlightCtrl);
// map.addControl(selectCtrl);
   highlightCtrl.activate();
// selectCtrl.activate();

then you get a working example, except that it does zoom to the cluster on hover.

I realy don't know how to fix this. Any ideas?
Many thanks in advance


texts: Where to get free seafloor elevation for Galapagos Islands?

texts:
I'm trying to find seafloor elevation around the Galapagos Islands. Something with resolution 30 - 100 meters.
Or just shaded relief for this region.


texts: Is table with 20,000 rows and 5 columns too large for Dissolve of ArcGIS for Desktop?

texts:
Is a table with 20 000 rows and 5 columns too large for arc map 10.2.1?

I have tried several things without success to dissolve this table based on one feature.

The table was obtained from the union of a shapefile table with a .xls table.

I can't even export the table to .dbf. When I do this, the resulting table only shows the first line.

Is this a size problem?


texts: PyQGIS method to add multiple GDB files?

texts:
I can use the following to add multiple shapefiles to QGIS:

import os, glob

layers=[]
shape_path = "C:\Users\Me\Desktop\Test\Shapefiles\\"
for shape in glob.glob(shape_path + "*.shp"):
    layer = QgsVectorLayer(shape, os.path.basename(shape), "ogr")
    layers.append(layer)
QgsMapLayerRegistry.instance().addMapLayers(layers)

However, I have dozens of .gdb folders which I would like to add automatically (I can do this individually through Layer > Add Vector Layer but would be quite tedious).

Is there a similar, or albeit another, method using PyQGIS to add multiple GDB files?


texts: I can't load a CSV file inside a QTableWidget by Python code

texts:
I am a beginner of PyQGIS and I'm trying to do a QGIS plugin, with QGIS 3.4.7 and PyCharm 2016.3.6. I already did the "dialog" box on Qt and after I wrote on Pycharm. It is a table only with the columns that I want.
I have the CSV file and my code. But when read the file directly on code and trying to write the iteraccion, it's kinda missing something, it means, the variable that belongs to the class QTableWidget (called tbwidget_data_list_service) on the code doesn't show up with .setItem (function and property of QTableWidget class). I have all libraries imported (PyQt5, and others), so I don't understand that's wrong.
I have some commented code 'cause I had many attempts to try this working...
Can someone tell me what's wrong or give another example?

Here's the code:
class DataListService(FrameworkSupportConsumer):

    def __init__(self, framework_support):
        super(DataListService, self).__init__()
        self.set_framework_support(framework_support)

        # Set QAction
        icon_path = ":/plugins/marineo_ocean/service_management.png"
        icon_action = QIcon(icon_path)
        message_action = 'Data List Service'
        self.action_data_list_service = QAction(icon_action, message_action, self.iface.mainWindow())

        # Get and register dialog
        self.dlg = DataListServiceDialog(self.iface.mainWindow())
        self.register_dialog('data list service', self.dlg)

        #  ler csv
        dadosaimportar = []
        ruta = "C:/Users/gviana/Desktop/MARINEO_PHASE_1/marineo_1/qgis/profiles/default/python/plugins/marineo_ocean/request.txt"
        reader = csv.reader(open(ruta), delimiter=';')
        for row in reader:
            dadosaimportar.append(row)
        numLinhas = len(dadosaimportar)
        numColunas = 7

        # self.tbwidget_data_list_service = QtWidgets.QTableWidget()

        #self.dlg.tbwidget_data_list_service.setRowCount(numLinhas)
        #self.dlg.tbwidget_data_list_service.setColumnCount(numColunas)

        for linha, umaLinha in enumerate(dadosaimportar):
            for coluna, value in enumerate(umaLinha):
                item = QtWidgets.QTableWidgetItem()
                item.setText(self,value)
                self.dlg.tbwidget_data_list_service.setItem(linha, coluna, item)

        # Connect events
        self.connect_signals()

    def unload(self):
        pass

    def connect_signals(self):
        self.action_data_list_service.triggered.connect(self.show_dialog)
        self.dlg.pushbtn_close_service.clicked.connect(self.close_dialog)

    def download_data(self):
        pass

    def display_data(self):
        pass

    def close_dialog(self):
        self.dlg.close()

    def show_dialog(self):
        self.dlg.show()


texts: Unicode problem in creating style with REST API on GeoServer

texts:
I am creating some new style with Postman. The main problem is that the new style was created with Unicode problem. It means I am seeing the new style has some character that are not Unicode.
URL : http://localhost:8090/geoserver/rest/workspaces/webgis/styles

HTTP Post

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
  <NamedLayer>
<Name>export_output8464daf0723</Name>
<UserStyle>
<Title>export_output8464daf071</Title>
<Name>export_output8464daf0732</Name>
 <FeatureTypeStyle>
   <Rule>
       <Name>export_output8464daf071_site_name_1</Name>
       <Title>export_output8464daf071_site_name_1_تپه قورچی</Title>
       <ogc:Filter>
         <ogc:PropertyIsEqualTo>
           <ogc:PropertyName>site_name</ogc:PropertyName>
           <ogc:Literal>تست</ogc:Literal>
         </ogc:PropertyIsEqualTo>
       </ogc:Filter>
       <PointSymbolizer>
          <Graphic>
         <Mark>
           <WellKnownName>circle</WellKnownName>
           <Fill>
             <CssParameter name="fill">#bd8321</CssParameter>
           </Fill>
         </Mark>
         <Size>2</Size>
       </Graphic>
     </PointSymbolizer>
   </Rule>
 </FeatureTypeStyle>
</UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

and Content type is :
application/vnd.ogc.sld+xml

Style in GeoServer is :
<?xml version="1.0" encoding="UTF-8"?><sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" version="1.0.0">
  <sld:NamedLayer>
    <sld:Name>export_output8464daf0732</sld:Name>
    <sld:UserStyle>
      <sld:Name>export_output8464daf0732</sld:Name>
      <sld:Title>export_output8464daf071</sld:Title>
      <sld:FeatureTypeStyle>
        <sld:Name>name</sld:Name>
        <sld:Rule>
          <sld:Name>export_output8464daf071_site_name_1</sld:Name>
          <sld:Title>export_output8464daf071_site_name_1_ط·ع¾ط¸آ¾ط¸â€، ط¸â€ڑط¸ث†ط·آ±ط¹â€ ط؛إ’</sld:Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>site_name</ogc:PropertyName>
              <ogc:Literal>ط·ع¾ط·آ³ط·ع¾</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <sld:PointSymbolizer>
            <sld:Graphic>
              <sld:Mark>
                <sld:WellKnownName>circle</sld:WellKnownName>
                <sld:Fill>
                  <sld:CssParameter name="fill">#bd8321</sld:CssParameter>
                </sld:Fill>
              </sld:Mark>
              <sld:Size>2</sld:Size>
            </sld:Graphic>
          </sld:PointSymbolizer>
        </sld:Rule>
      </sld:FeatureTypeStyle>
    </sld:UserStyle>
  </sld:NamedLayer>
</sld:StyledLayerDescriptor>


texts: How to convert masks generated by rasterio to channel last mask

texts:
The problem is about using rasterio to generate mask from GeoJSON file which contains geographic coordinates.

The function that I used is from:

from rasterio.tools.mask import mask
with rasterio.open(img_path) as src:
        img_mask, out_transform = mask(src, [poly], crop = False)

The shape of the mask that I generated is:

(4, 256, 256)

I converted this mask to (256,256) black and white image by taking the first channel of the mask. I used for loop to do this process couple times, and then I appended all the (256,256) masks into a list. I tried to use np.stack to convert the list to a numpy array with dimension (256,256,?). However, I got an array with dimension (?,256,256). I tried to set axis=-1; however, the image shape becomes (256,?) rather than (256,256).

Is there any method that can convert the channel first array to a channel last array with numpy?