texts:
I'm trying to create a cluster with categories using PruneCluster around my GeoJSON layers.
As the example in GitHub They are creating randomly some markers then adding them to the cluster.
I did the same, I got the coords from a GeoJSON layer, stored it in a variable then added it to Prunecluster as follows :
leafletView.BuildLeafletClusterIcon = function(cluster) {
        var e = new L.Icon.MarkerCluster();

        e.stats = cluster.stats;
        e.population = cluster.population;
        return e;
    };

    var colors = ['rgb(10, 94, 219)', 'rgb(253, 152, 1)', 'rgb(3, 180, 11)'],
        pi2 = Math.PI * 2;

    L.Icon.MarkerCluster = L.Icon.extend({
        options: {
            iconSize: new L.Point(44, 44),
            className: 'prunecluster leaflet-markercluster-icon'
        },

        createIcon: function() {
            // based on L.Icon.Canvas from shramov/leaflet-plugins (BSD licence)
            var e = document.createElement('canvas');
            this._setIconStyles(e, 'icon');
            var s = this.options.iconSize;
            e.width = s.x;
            e.height = s.y;
            this.draw(e.getContext('2d'), s.x, s.y);
            return e;
        },

        createShadow: function() {
            return null;
        },

        draw: function(canvas, width, height) {

            var lol = 0;

            var start = 0;
            for (var i = 0, l = colors.length; i < l; ++i) {

                var size = this.stats[i] / this.population;

                if (size > 0) {
                    canvas.beginPath();
                    canvas.moveTo(22, 22);
                    canvas.fillStyle = colors[i];
                    var from = start + 0.14,
                        to = start + size * pi2;

                    if (to < from) {
                        from = start;
                    }
                    canvas.arc(22, 22, 22, from, to);

                    start = start + size * pi2;
                    canvas.lineTo(22, 22);
                    canvas.fill();
                    canvas.closePath();
                }

            }

            canvas.beginPath();
            canvas.fillStyle = 'white';
            canvas.arc(22, 22, 18, 0, Math.PI * 2);
            canvas.fill();
            canvas.closePath();

            canvas.fillStyle = '#555';
            canvas.textAlign = 'center';
            canvas.textBaseline = 'middle';
            canvas.font = 'bold 12px sans-serif';

            canvas.fillText(this.population, 22, 22, 40);
        }
    });

let myLayerOptionsArbre = {

        pointToLayer: createCustomIcon,
        onEachFeature: function(feature, layer) {
            layer.bindTooltip(feature.properties.titre);
            layer.on({
                click: onMarkerClick
            });
            coordsArb.push(feature.geometry.coordinates);

        }
    }

var arb = arbres
    var arbre = new L.GeoJSON(arb, myLayerOptionsArbre)

var ap = {}

    for (var i = 0; i < coordsAP.length; i++) {

        ap = new PruneCluster.Marker(coordsAP[i][1], coordsAP[i][0]);
        console.log(typeof ap)
        console.log("hey", coordsAP[i][1])
        ap.category = AIRES_PROTEGEES;
        leafletView.RegisterMarker(ap);
    }

leafletView.ProcessView();
    mymap.addLayer(leafletView)

Now on my map I have the clusters around new markers and the GeoJSON layers. When I read the doc, I found that Prunecluster does not support GeoJSON.
Now if I just work with markers, How can I attach them to the properties of the geojson features, because I'm building custom popups with the function onMarkerClick called when a geojson feature is clicked.
So if you have an idea about how to achieve that I would be pleased.


texts: Making visual realistic map containing mountain hiking trails/tracks?

texts:
I would like to try making a map that visually reproduces a mountain/hill and contains hiking trails, just like this one:
http://i.stack.imgur.com/PdcnM.jpg

I imagine some steps like this:

Generate a realistic 45 degrees (or so) representation of the area (from GIS) - how?
An artist (painter) draws the visual simplified representation of what you get on firsts step and then trails are added.

Did anyone do maps like this and knows exactly how to make one?


texts: What CRS to use when creating buffer all over the world in QGIS

texts:
I trying to create buffers (50km) around many points all over the world (airport locations) in QGIS. I will then overlay these buffers with a world-wide raster dataset so I can extract information from that. The point layer is in WGS84 now, I am wondering what is an appropriate projected coordinate to use (one that does not give me too much distortion)?

I have searched this site and the internet a lot and did not find an answer because I think this cannot achieved with PostGIS or just by using a projected coordinate system for one country.


texts:
What's the best way to determine if a point is within a certain distance of a GEOJSON polygon? Should one use TurfJS buffer method (https://github.com/Turfjs/turf-buffer#turf-buffer)? Can one perform queries on the buffered polygon?

It's clear to me one could us the TurfJS' inside method (https://github.com/Turfjs/turf-inside) to determine whether a point is within a polygon. I'm just curious what the best approach would be for finding whether or not a point is inside of a buffered polygon.

For example: I have a number of neighborhoods provided as a GEOJSON polygon files. I also have a set of locations/addresses for employees (already geocoded to lat/long coordinates). What would be the best way to see whether or not my employees live within 10 miles of a given neighborhood polygon?


texts: Changing Layer DataSource using ArcObjects?

texts:
I want to change data source information of a layer. Lets say I have a MXD file in dev env configured with dev environment datasource details like server,database,username,password and I want to use this MXD in test environment pointing to test datasource details. I have created a Add-In button. On click on this button, ArcMap Database Connection windows pops up asking database details to which layer should point in test environment. Once user add the details I capture this information in IWorkspace workspace:

//Code for popping up a window to ask user for datasource details

IPropertySet propertySet = new PropertySetClass();
propertySet.SetProperty("SERVERINSTANCE", "sde:sqlserver:devserver"); //sde:sqlserver:(servername)
propertySet.SetProperty("dbclient", "SQLServer");
propertySet.SetProperty("DATABASE", "databasename"); //(databasename)
Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SdeWorkspaceFactory");
IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);
IWorkspace targetWorkspace = workspaceFactory.Open(propertySet, 0);

In below code, I am trying to access layer's workspace connection properties and then trying to set it. I dont know how to proceed further with the code.

//Code for trying to set layer datasource.

IDocument doc = ArcMap.Application.Document;
IMxDocument mxDoc = doc as IMxDocument;
IMap map = mxDoc.FocusMap;
ILayer layer = map.get_Layer(0);// Layer[0];
ESRI.ArcGIS.Carto.IFeatureLayer featureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)layer;
ESRI.ArcGIS.Geodatabase.IFeatureClass featureClass = featureLayer.FeatureClass;
ESRI.ArcGIS.Geodatabase.IDataset dataSet = (ESRI.ArcGIS.Geodatabase.IDataset)featureClass;
ESRI.ArcGIS.esriSystem.IPropertySet propertySet = dataSet.Workspace.ConnectionProperties;

dataSet.Workspace.ConnectionProperties.SetProperty("Server", "testservername");


texts: I am looking for Thermal IR Photography with 1 m resolution?

texts:
I am looking for thermal IR photography with 1m resolution?

It needs to be recent or maybe a company in western Canada that can fly a mission like this. Satellite or flown photos both work, I am just looking for sources and maybe the specific satellites to look for.


texts:
I have a huge area of linestrings and I want to get the features of each portion of area. All the linestring parts are stored in a python list and they are not in order.

One portion of area is shown in the figure 1.

I want to create a bounded box on specific percentage of my area and then get the feature out of each linestring which is intersecting with bounded box. This is illustrated in figure 2. The red box is bounded box on the linestring and it should get the percentage of each intersecting line

My problem is that I also want it to slide to the whole area from one starting point. I have seen that python shapely has a box but dont know how to slide it to the whole area. Figure 2 and 3 shows how I want sliding..

This is simple illustration in reality I want a box that slides from start of my area to the end of my area which is stored in python list.

Shows bounded box is sliding on the area

Figure 1 to 3 show how it should be sliding.

Is it possible to implement in python? I want to get features of both side of linestring which is a road or is there some other method?


texts: Converting Landsat8 DN to TOA reflectance for use in envi

texts:
I am used to working in envi, however my current version does not support Landsat8.
I downloaded GRASS for the express purpose of converting Landsat8 DN to TOA reflectance for use/analysis in envi.

I have managed to do the conversion fine using the tutorial provided (http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/#Unpacking_the_downloaded_Landsat_8_dataset), however once I try to export that to a format usable by envi, I get stuck.

I have tried "r.out.tiff":
However get this warning:
"WARNING: Raster map > is a floating point map. Fractional values will be rounded to integer"

When I load it in envi it is clear that the image quality has decreased and it is no longer suitable for my analyses.

When I try "r.out.gdal":
I get the following error message: "ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format."

And when I try to open it in envi, I get this message: "Cannot open in envi: unsupported number of bits per pixel"

Can anyone help me please?


texts: Using CartoDB Infowindow on fullscreen of iPad?

texts:
We're just about ready to go with our first ever CartoDB visualisation, namely a map of the museums in Britain.

Unfortunately I'm having a rather large issue with infowindow positioning in either desktop fullscreen or when viewed on an iPad. In short, when clicked on fullscreen then the infowindows appear about 1,000 miles south of their marker. On an iPad it seems to be around 1,000 miles north of their marker. The code for the infowindow is as follows:

x

<divclass="cartodb-popup-content">

{{destination_name}}

county

{{county}}

Type of Museum

{{destination_type}}

description

{{description}}

Visit the Museum's Website

The URL for the visualisation is at http://www.historic-uk.com/HistoryMagazine/DestinationsUK/Museums/


texts:
In QGIS WMS/WMST connection we can define a Referer (URL) as can be seen here:

Since I was not able to find an exact definition for this:

Q: What is a referer and what is it used for?


texts: How can I represent distance using the dashes or markers of a vector line?

texts:
I am a very basic user and dont know the lingo to get to the gist of my question so some background might help.

Ive volunteered to create a pre-incident planning system for my rural fire department. Our primary vulnerability are extremely long un-mapped driveways.

So far Ive successfully transferred GPS data of the driveways over aerial photography of the district as bright green lines(a .SHP layer) with a label showing the length of each driveway. It looks great and is very useful as is.

Now the question: I want to draw the driveway vector lines as dashes or with markers set at 100-foot lengths. Our supply-line hoses come in 100-foot lengths; A vector line defined with markers at 100-foot intervals would save me a lot of time and improve the data immeasurably.

So far I have tried "Single Symbol" as a "Simple Line" using a customized dash pattern and the "Marker Line" using the interval placement. Neither of these options maintain the proper length while zooming or "scaling" (if that's the term).

Perhaps it is impossible? (at least at my level of knowledge. I taught myself this far and dont have a lot more time for fooling around :))


texts: Setting top-level GeoJSON id using ogr2ogr?

texts:
I'd like to write the contents of a PostGIS table (or query) to GeoJSON using ogr2ogr.

How can I select an attribute and have it be used as a top-level id in the output GeoJSON, without losing it as a property?

Given some PostGIS data:

CREATE TABLE things (
    city_id int,
    entity_id int,
    geom geometry(Point, 4326)
);

INSERT INTO things VALUES (3, 19, 'SRID=4326;POINT (4 8)');

And an ogr2ogr command:

ogr2ogr -f "GeoJSON" things.geojson PG:"dbname=postgres" ***magic*** things

I'd like to see:

{
    "type": "FeatureCollection",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
        }
    },
    "features": [
    {
        "type": "Feature",
        "id" : 19,
        "properties": {
            "city_id": 3,
            "entity_id": 19
        },
        "geometry": {
             "type": "Point",
              "coordinates": [ 4.0, 8.0 ]
        }
   }]
}

Note that the value of entity_id was copied to the top level, but the attribute was not removed from the properties section.

I've tried the following:

Set primary key to entity_id then run ogr2ogr with -preserve_fid.  This gives me a top-level "id" but removes entity_id from the properties.
Use SQL to copy the entity_id to another column, and then instruct ogr2ogr to use this as the id.

ogr2ogr -f "GeoJSON" things.geojson PG:"dbname=postgres"
--config PGSQL_OGR_FID=fizz
-preserve_fid
-sql "SELECT *, entity_id AS fizz FROM things"

In this case, I get no top-level id.

Tried this on GDAL 1.11.2 and current master.


texts: Add .OSM (OpenStreetMap) data to a QGIS project using python

texts:
I try to load a .osm file into my QGIS project using python.
So I've write this code:
import processing
from qgis.core import *

# Creation of the .osm file, the QuickOSM plugin is needed.
bbox = '795078.7779580103, 795572.1303964029, 5908035.195114954, 5908380.196820124'

alg_params = {
    'EXTENT': bbox,
    'KEY': 'highway',
    'SERVER': 'https://lz4.overpass-api.de/api/interpreter',
    'TIMEOUT': 25,
    'VALUE': ''
}

query = processing.run('quickosm:buildqueryextent', alg_params)
file = processing.run("native:filedownloader", {'URL':query['OUTPUT_URL'], 'OUTPUT':'TEMPORARY_OUTPUT'})

# Now I try to display the osm data:
vlayer = iface.addVectorLayer(file['OUTPUT'], "lines", "ogr")

The code do not crash, but since an .osm file can contains several different geometry the iface.addVectorLayer() function open this windows:

And I need to manually select the geometry that I want to display (lines, multipolygons...). Is there a way to bypass this windows ? Is there a way to either load all the different geometry or to select the geometry that I'm interested in directly in the iface.addVectorLayer() function ?


texts: Geoexplorer displays the layers I import, but not the table behind them

texts: If I import layers into GeoSrver, they display in GeoExplorer, and I can use the identify tool to find more info about areas of the layer I select. However the table from which this data comes from does not show at the bottom of the screen, as it does for existing example layers. I assume there is a setting somewhere to enable this?


texts: Colour issue when reprojecting image mosaic

texts:
We have a couple of instances of GeoServer running, one of which is the latest version and the other is version 2.1.3.

The older version has Ordnance Survey base mapping which is loaded as a number of EPSG:27700 image-mosaics. These are then reprojected by GeoServer into the Google projection EPSG:900913 and displayed as part of a layer group where the style switches the layers on / off depending on the scale. We believe this is a fairly standard use of GeoServer to provide a base map coverage and it works fine.

In the newest versions of GeoServer the same base map layer creates a horrible orange coloured map as it seems to be breaking the colour banding of the layers. Seems to be turning the large swaths of blue sea into orange, as well as screwing the colours on the rest of the map. I can't upload a screenshot as our organisations internet policies restrict it.

Please note, you need to reproject two image-mosaics and combine them in a layer group or GetMap to reproduce the error.

We wondered if anyone else has experienced the bug in the latest versions, and whether there are any solutions or workarounds anyone has come across.


texts: No response from addPart operation when trying to upload parts of .tpk to AGOL (Python, requests, ArcGIS REST API)

texts:
Trying to automate the updating of a large .tpk on AGOL. When the .tpk is updated on my servers I do the following:

Break into parts
Call updateItem operation with multipart = true
Make a request to upload each part using addPart
Make a request with the commit function to finalize the update to the .tpk on AGOL

My issue is that I get no response when trying to use addPart, then when the commit operation is hit in my code it says there are no parts to commit. A very odd result as something is definitely happening when the addParts operation is performed. Memory usage kicks up to ~ 3GB (1gb files so this seems right) and each addPart takes maybe 5 minutes or so complete and move to the next. It's not a firewall/proxy issue with the request as I've made multiple request in this environment to AGOL and had no issues getting responses.

Code:

update_url = 'https://lipa.maps.arcgis.com/sharing/rest/content/users/{}/items/{}/update'.format(username,tpk_id)

up_data = {'f':'json',
    'token':token,
    'title': layer_name,
    'itemId':tpk_id,
    'type': 'Tile Package',
    'overwrite': 'true', # existing item, overwrite is set to true
    'async':'true',
    'multipart':"true"}

update_response = requests.post(update_url, up_data=data, verify = cert).json()

print update_response

part_num = 1
part_url = 'https://lipa.maps.arcgis.com/sharing/rest/content/users/{}/items/{}/addPart'.format(username,tpk_id)

for part in os.listdir(parts):

    path = r"\\psegliny.com\oms\oms_gis_prd\fileshare_mapping\cache_extent\es_secondary_tpks\full_service_area"

    part = path + "\{}".format(part)

    print part
    print(str(part_num))

    ##files_up = {"file": open(part, 'rb')}
    with open(part, "rb") as f:
        files_up = {"file":f}

        data = {'f':'json',
                'token':token,
                'partNum':part_num}

        part_response = requests.post(part_url,data = data,files=files_up,verify=cert)
    part_num+=1
    try:
        jres = json.loads(part_response.text)
        print jres
    except:
        print "no response"
    ##print add_part_response

commit_url = 'https://lipa.maps.arcgis.com/sharing/rest/content/users/{}/items/{}/commit'.format(username,tpk_id)

data = {'f':'json',
'token':token}

commit_response = requests.post(commit_url,data=data,verify=cert).json()

print commit_response

Anyone know what I am doing wrong? It's very odd that I get no response at all, at the very least I should get a failed/success response back (as per the documentation)


texts: Deleting edge pixels in raster using ArcMap

texts:
Does anyone know how to delete single edge pixels in a raster after a mask is completed?

I could convert to polygon and back to raster, but I don't want to lose accuracy. I am using ArcMap.


texts: Add delimited text layer is no longer an option in Data Source Manager for QGIS 3.10.6

texts: I have a CSV with lat/long coordinates that I wish to turn into a shapefile. The process involves going to Layer > Add Layer > Add Delimited Text Layer. In QGIS 3.4.6, this would open up the Data Source Manager, where on the right-hand column there was a "Delimited Text" option to import a CSV and create a shapefile (First image). I recently installed QGIS 3.10.6, opening Data Source Manager no longer includes a "Delimited Text" option to import a CSV with lat/long coordinates (Second image). Is there a way to add this option to the Data Source Manager? Or do I have to find another way to accomplish this task.


texts:
Newish to R
Very new to GIS

I've grabbed some vector spatial census data using the cancensus package.
I now have a spatial data frame with variables and geometry
I plot it and decide to focus on a part of the whole province

My thinking here is to create a rectangle with the coordinates that I am interested in using extent(xmin, xmax, ymin, ymax)

Assign it the same CRS as the spatial data frame and then use st_intersection. Is this the right approach?

extent(box) <- c(-81, -74, 41.68, 45)
st_set_crs(box, "+proj=longlat +datum=WGS84 +no_defs") # Fails
ont_crop <- st_intersection(census_data, box)

I'm not sure how to show you a reproducible example since census data is large. Please let me know what I am doing wrong.


texts: LDAP user not accessible published layers

texts:
I’m facing an issue regarding the roles associated with the LDAP users in Windows Server 2019
Below steps have been followed

Created the LDAP group and user with organizational unit associated them with defined user  in Windows Server active directory.

I have created authentication filter in GeoServer and the connection was successful with LDAP user using GeoServer and was able to login it in GeoServer

I have created successfully roles in user group service and selected the LDAP defined role in GeoServer

Now I have added few layers and associated them with LDAP user in GeoServer.

The issue I’m facing those layers I’m not able view with LDAP viewer whereas I can view with admin user all the layers
Where did I make the error?


texts: Tag for calculating way length in osm2pgsql?

texts: There is a special tag 'way_area' in osm2pgsql creating a database column that stores the area of polygon features. So is there something like this for calculation of length of line objects?


texts: How to keep selected layers on page load in Leaflet?

texts:
I have a leaflet map that reload page and map view depending on a drop-down list. I have many layers in it and I'd like to keep selected layers as selected when the page reload after clicking on an item in the drop-down list.

I hope it's clear enough.


texts: How to make 2 raster to be the same extent?

texts:
I got two raster files of rice suitability area and land use. Two layers have the same cell size (400x400 meters). However, rice suitability layer extent is (234,298) column and row, while extent of land use layer is (229,288). It caused a big trouble when I run allocation model. I would like to make extent of suitability layer to be the same with land use layer.

I have tried to use data management> raster> raster processing> clip and selected environment> processing extent> select extent same as layer land use and/or snap raster to "land use" layer. Unfortunately, the extent of rice suitability layer is the same.

How can I fix this problem with ArcGIS?


texts: How to resolve this error 010067 while running a reclassified Raster to Polyline?

texts: I'm facing this error while running a reclassified Raster to Polyline in ArcGIS 10.2. Also, I've attached the input data information


texts:
I have a leaflet map in a bootstrap 3 html application. I have content both above and below the the map which has a height of 600px. When I access the application from a phone, the width of the map adjusts to the the width of the screen. However, the height of the map does not update and you end up with a seemingly endless map sandwiched between two chunks of content. The user must then endure endless scrolling. You can see what i mean here: http://punkmap.com/examples/fauxnmap/index.html

Anyone know what's the best way to handle this?


texts: Changing default decimal places displayed in ArcGIS Pro

texts:
I am using ArcGIS Pro 1.4.0. My current problem is that all new numeric values (float and double) are being rounded into whole numbers in the attribute table of shapefiles. According to this ArcGIS Desktop help site Formatting numeric fields in tables ,

  "six is the default number of decimal places to display in the Table window."

However, this does not seem to be the case for my version of ArcGIS. For example, see below image. First, I calculated the geometry of gridcells in my shapefile (POLY_AREA). Next, I created a new field and calculated it to equal POLY_AREA. As displayed, all numbers are rounded to the whole number. I've included the python window in the clipped image for reference. Also, I am not sure what version of ArcMAP the original shapefiles were created in, but definitely pre-ArcGIS pro. The issue continues to occur even after exporting/copying them several times within pro.

In this instance, the rounding is not a problem; however, moving forward I will be conducting calculations on many small decimals. I know how to manually change the display, but it is very tedious to do tens or hundreds of times.

Is this a bug, or does anyone know how to globally change the default display of numeric values?


texts: Create Map Tile Package (tpk) progress? arcpy.CreateMapTilePackage_management

texts:
I'm using arcpy.CreateMapTilePackage_management to happily create a .tpk - all good.

Trouble is this can take a matter of hours, or days for a large complex mxd, once arcpy.CreateMapTilePackage_management has been called, is there any way of getting progress of how it's doing?


texts: QGIS 3 Value Relation current_value() not working in custom ui form

texts:
I have a custom UI form made with Qt Designer with a couple of dropdown list, so I can easily edit a layer. I use a Value Relation on those fields to filter the displayed data. The problem is the function current_value() is not working when I use a provided ui file, but it works perfectly when I use the default autogenerated form. I've tried this both programatically and through the Layer Properties interface with the same result. I guess the problem is the function current_value() doesn't find the field to be filtered by, but I don't know how to point to that field correctly on my custom form. I want to set the state first, and then show only the towns that are on that state.

The form widget name is the same as the field name. I tried this both in QGIS 3.4.3 and 3.6.1. This is the code I used, though it worked fine on default form. Im I missing something or is this a bug ?

myLayer = QgsVectorLayer(uri.uri(), 'myLayer', 'postgres')
states = QgsVectorLayer(uri.uri(), 'states', 'postgres')
towns = QgsVectorLayer(uri.uri(), 'towns', 'postgres')

stateIndex = myLayer.fields().indexFromName('state')
townIndex = myLayer.fields().indexFromName('town')

myLayer.setEditorWidgetSetup(stateIndex, QgsEditorWidgetSetup('ValueRelation', {
    'Layer': states.id(),
    'Key': 'id',
    'Value': 'name'
    }
))

myLayer.setEditorWidgetSetup(townIndex, QgsEditorWidgetSetup('ValueRelation', {
    'Layer': towns.id(),
    'Key': 'id',
    'Value': 'name',
    'FilterExpression': '"fk_state" = current_value(\'state\')'
    }
))


texts: Get source and output paths in QGIS 3 processing algorithm

texts:
How can I find the complete path of the INPUT (r_in_name) and OUTPUT (r_out_name) rasters for the gdal datasets in the following code:
from qgis.processing import alg
from qgis.core import QgsFeature, QgsFeatureSink
import numpy as np
from osgeo import gdal, gdal_array, osr
@alg(name="quantilizearaster", label=alg.tr("Quantilize a raster)"), group="rastertools", group_label=alg.tr("Raster Tools"))
@alg.input(type=alg.RASTER_LAYER, name="INPUT", label="Input layer")
@alg.input(type=alg.RASTER_LAYER, name="OUTPUT", label="Output layer")
def quantilize(instance, parameters, context, feedback, inputs):
    """
    This algorithm takes an INPUT raster and quantilize its values in five classes; result is a quantilized OUTPUT raster.
    """

source = instance.parameterAsRasterLayer(parameters, "INPUT", context)
output_raster = instance.parameterAsRasterLayer(parameters, "OUTPUT", context)

r_in_name =   # r_in_name must be a string of the complete path to the INPUT raster, including its file name (*.tif)
r_out_name =  # r_out_name may be a temporary OUTPUT raster file

dataset = gdal.Open(r_in_name)
band = dataset.GetRasterBand(1)
nodata = band.GetNoDataValue()
array = dataset.ReadAsArray()

new_array = array
nan_array = array

nan_array[array == nodata] = np.nan

percentile_80 = np.nanpercentile(nan_array, 80)
percentile_60 = np.nanpercentile(nan_array, 60)
percentile_40 = np.nanpercentile(nan_array, 40)
percentile_20 = np.nanpercentile(nan_array, 20)
percentile_0 = np.nanpercentile(nan_array, 0)

for i, v in enumerate(new_array):
    for j, element in enumerate(v):
        if element > percentile_0 and element <= percentile_20:
            new_array[i,j] = 1
        if element > percentile_20 and element <= percentile_40:
            new_array[i,j] = 2
        if element > percentile_40 and element <= percentile_60:
            new_array[i,j] = 3
        if element > percentile_60 and element <= percentile_80:
            new_array[i,j] = 4
        if element > percentile_80:
            new_array[i,j] = 5

new_array[ new_array != new_array ] = nodata

geotransform = dataset.GetGeoTransform()
wkt = dataset.GetProjection()

# Create gtif file

driver = gdal.GetDriverByName("GTiff")
output_file = r_out_name  # r_out_name may be a temporary OUTPUT raster file

dst_ds = driver.Create(output_file,
                       band.XSize,
                       band.YSize,
                       1,
                       gdal.GDT_Int16)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( new_array )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(nodata)
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset

dataset = None
dst_ds = None

return


texts: How to create a chart of elevation using points?

texts:
I am using shapefiles with basic information including the x-, y- and z-coordinates. I want to see how the objects (in this case points, which are finds from an excavation) are scattered along the z-axis. To make it visible, I copied the information to Excel and made a chart. There I can see the distribution of the finds on the z-axis, of course. But I want to make this distribution visible in QGIS (with the z-axis and x/y-axis) like as it can be seen in charts made with Excel.

I already tried to do this with the tools "Profile from line" and "Profile tool", but I always failed because I neither use raster files nor are the contents of the shapefiles drawn by lines, they are only drawn by points.

Is this possible with QGIS? If so, how is it done? (I'm using QGIS Lisboa)


texts: Downloading GeoExplorer and installing outside suite?

texts:
I already have QGIS, GeoServer and QGIS all configured, but I would love to add GeoExplorer so I can publish maps easily.

Is there a way to do this without having to install the whole opengeo suite?

Is it possible to install GeoExplorer on its own to work with QGIS?

Are there other options for publishing as a web application with edit functionality?


texts: TRMM Values Conversion From Float to Integer in Google Earth Engine?

texts:
In the below code TRMM product values converted from float to integer type in Google Earth Engine. However after the conversion TRMM product just shows 0 and 1. How can I can do conversion without this error?
code link: https://code.earthengine.google.com/76c3540de4e0fe14d022dc1bf789d586
Map.centerObject(table);
Map.addLayer(table);

var start = '2018-01-01';
var end = '2019-01-01';

var precipitation = ee.ImageCollection("NOAA/PERSIANN-CDR")
.filterDate(start, end)
.filterBounds(table)
.select('precipitation').mean().toInt16()
.clip(table)
.rename('precipitation');

print('precipitation',precipitation);

Map.addLayer(precipitation,{
  palette: ['blue', 'purple', 'cyan', 'green', 'yellow', 'red']
},'precipitation',false);


texts: WFS GetFeature request with CQL_FILTERS with SPATIAL OPERATORS does not retrieve correct results

texts:
I have a WFS service :

https://services.ovam.be/geoserver/wfs?REQUEST=GetCapabilities&SERVICE=WFS

Where I'm trying the following:

http://services.ovam.be/geoserver/wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=9999&CRS=EPSG:31370&BBOX=123604.6385,210458.6915,123764.6385,210618.6915&OUTPUTFORMAT=json&TYPENAME=OVAM:Dossierinfo&PROPERTYNAME=(gid,dossier_id,dossier_obo,dossier_bbo,dossier_bsp,dossier_eeo,dossier_status,status_jaar,dossier_type)

This gets me the results within the BBOX.
This includes Features in the whole BBOX.
I only want to get the Features within a specific range.

My spatial box (the container I want results for) is a Polygon :

<?xml version="1.0" encoding="UTF-8"?>
<Polygon xmlns="http://www.opengis.net/gml" srsName="http://www.opengis.net/def/crs/EPSG/0/31370">
   <exterior>
      <LinearRing>
         <posList>123664.03100000322 210489.5150000006 123705.24599999934 210489.5150000006 123705.24599999934 210587.86800000072 123664.03100000322 210587.86800000072 123664.03100000322 210489.5150000006</posList>
      </LinearRing>
   </exterior>
</Polygon>

I have tried building a query using the CQL Filter:

CQL_FILTER=INTERSECTS(geom, POLYGON((123664.03100000322 210489.5150000006,123705.24599999934 210489.5150000006,123705.24599999934 210587.86800000072,123664.03100000322 210587.86800000072,123664.03100000322 210489.5150000006)))

Full URL :

http://services.ovam.be/geoserver/wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=9999&CRS=EPSG:31370&BBOXs=123604.6385,210458.6915,123764.6385,210618.6915&OUTPUTFORMAT=json&TYPENAME=OVAM:Dossierinfo&PROPERTYNAME=(gid,dossier_id,dossier_obo,dossier_bbo,dossier_bsp,dossier_eeo,dossier_status,status_jaar,dossier_type)&CQL_FILTER=INTERSECTS(geom, POLYGON((123664.03100000322 210489.5150000006,123705.24599999934 210489.5150000006,123705.24599999934 210587.86800000072,123664.03100000322 210587.86800000072,123664.03100000322 210489.5150000006)))

The resultset is empty - although there should be results ...
If I call the GetFeature request having a BBOX of 123664.03100000322,210489.5150000006,123705.24599999934,210587.86800000072 it does retrieve results ... This BBOX is the lower / highest spatial point in the above POLYGON definition.

http://services.ovam.be/geoserver/wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=9999&BBOX=123664.03100000322,210489.5150000006,123705.24599999934,210587.86800000072&OUTPUTFORMAT=json&TYPENAME=OVAM:Dossierinfo&PROPERTYNAME=(gid,dossier_id,dossier_obo,dossier_bbo,dossier_bsp,dossier_eeo,dossier_status,status_jaar,dossier_type)&CRS=EPSG:31370

Any help in what I'm doing wrong?


texts: Restricting loaded area in Autodesk Map3d for WMS connection?

texts:
I have an in house WMS server running on Ubuntu Linux and using MapServer to serve a WMS layer for my various bits of GIS software.

The server works perfectly fine, but it does take a little while to create a map if it's very large.

I have a tile layer for example that covers a very large area of the north east of England, from Blyth right down to North Yorks, and from the Pennines out past the East Coast.

Loading this WMS layer into QGis (as an example) is perfectly fine, as long as I Zoom in, and don't zoom out too far, and since I'm using this as background mapping to see where I need to draw things, then it's not an issue.  Mapserver produces the map, for the area requested.

My problem with all of this, is when I use Autodesk Map3D.

It accepts the WMS connection fine, but no matter how far in I'm zoomed, every time I try to add the WMS connection it tries to add the entire map into the Map3D project.  This (because of the size of the map) times out every time, thus preventing me from adding the background map to my project.

Does anyone here know how I can get Map3d to only load the area under the view port, then load the other parts as needed when I scroll.

As I say, QGis seems to do this perfectly fine, and even my Open Layers web viewer Iv'e written works perfectly too, but Map3d just will not, and right now I can't seen to find any sensible answers anywhere.


texts:
I'm loading gtfs into a PostGIS database but I have a problem when loading my stop times.

The data is in time format "00:00:00" and for 1:30 AM, my data is "25:30:00".

I'd like to keep it that way because usually, night buses services start before midnight so I can keep track of the schedule by time, but PostgreSQL doesn't accept time after 24:00:00.

I could modify my script (python) so when it comes to treating stop_times.txt, to modify the time to create something like 2530 so I'll have to deal with integer, but I'm not the only one  that will have to deal with that db and my colleague would really prefer to deal with time instead of integers...

Does anyone could think of a way to "bypass" this limitation?

Thanks!!


texts: How to convert web-map website to mobile phone

texts:
I created simple web-map application,It works perfect on Desktop,but when I open on mobile or tablet then their is size problem.

When I open web-map in mobile then all the web-map content size are changed. I also added this html code

<meta name="viewport" content="width=device-width, initial-scale=1.0,maximum-scale=1.0,user-scalable=no">

in my file but still not showing proper output.

How I convert web-map to responsive web-map.How to test it on local machine.


texts: How to control which Fields are exported in ArcObjects? (via Field Mapping)

texts:
In ArcMap 10.x, if you use the Export Data command on a layer that has hidden/disabled fields, those fields will not be present in the created feature class/table. My question is simple: how do you duplicate this functionality using ArcObjects and C#? The IField interface seems to have no property or method for defining this aspect of table data. Can it be done?

I am seeking to do this so I can extract a subset of fields from a table and keep it for later use (as opposed to having to deal with all the unnecessary data bloat present in many standardized data releases).

UPDATE: I am currently attempting to use the IFieldInfo option suggested in the comments below, but I find myself unable to perform the conversion from an IFeatureClass to an IFeatureLayer that seems necessary to access ITableFields. Am I missing something? I am doing this in a WinForms app, so there is no parent ArcMap instance open to obtain layer information from. I wish to only work with feature classes and tables, if possible.

UPDATE 2: It would seem that, because this winforms is detached from any instance of ArcMap, I cannot make use of IFeatureLayer to turn fields on and off before exporting. Does anyone know how to properly create a field map object in C#?


texts: Define custom projection with small offset in QGIS

texts:
I have a series of raster maps of geophysics for an archaeological site in Turkey. These are ostensibly in UTM zone 35, WGS84 projection, but the surveyors set their grid with a total station using a slightly different coordinate for our local datum than the one we use, 1.726 m to the W and 0.735 m to the S. So when I bring them into QGIS 3.18 with the data for the site as a whole, they end up shifted a couple meters to the SW. What is the easiest way to bring them all into alignment with the rest of the system?
I naively tried creating a custom projection with a small offset (like post Shift a CRS by an observed offset) by adding 1.726 to the false easting and 0.735 to the false northing of the project CRS 32635, so:
PROJCRS["GeoIm 2018",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",27,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",50001.726,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0.735,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
    ID["EPSG",9001]]]]

This however shifts a test raster by about 450000 m to the east. I'm obviously flailing here; can someone point me in the right direction? Since I will be getting more data from these folks in the future, it seems better to somehow recognize their system than manually move every raster by a given amount.


texts: I have a NDVI image, which I Geo-referenced in QGIS (raster).  Now I have to open the image in Google Earth - how do I convert the GeoTIFF raster file to a kml file in QGIS, so that I can open the image as a layer in Google Earth?


texts: Does CARTO Airship provide a feature for widget attribute editing?

texts:
CARTO Builder allows users to edit feature attributes directly from the map building and editing interface.

Does Airship offer the same function?


texts:
I have a polyline feature class and a polygon feature class. I would like to clip the extent of the polyline which are NOT overlaid on the polygon.

How can I do that?

I should be looking for similar tool as intersect, but intersect clips the polyline that overlaid on the polygon. I want an opposite of that.


texts:
I'm working with QGIS Brighton. My project was started with an older version of the software, but I don't think that makes a difference. I have six points that were generated  by calculating the intersection of other layers that are no longer part of my project. At the time, I added a pinned label to each point. I thought that when I would delete the layers that created the intersections that the labeled points would also disappear. This didn't happen and I can't seem to delete the pins no matter what I do since they don't exist on a layer that is in the project.

I'm attaching a picture of the pinned point (circled in blue).


texts:
I'm new to GIS and mapping services.

Can anyone tell me what my options are for reverse-geocoding libraries, other than nominatim? I'm most likely going to be working with OpenStreetMap data. So I'll be on the PGSQL/PostGIS stack. Does any of them support the notion of "radius" in their search (which standard Nominatim doesn't support)?


texts: {z}/{x}/{y} tiles for Belgium names only

texts:
Right now I have a leaflet map, with GEOJSON layers rendered on top of them, behind that I render {z}/{x}/{y} tiles that only displays names of countries, cities, and so on.
The thing is, that everything that falls outside Belgium, i don't really want to see on the map.
I searched a lot but can't find one with only names

That why I wondering if there are any free name layers from ONLY belgium on the internet?
If not, is it hard to make one?


texts: Rendering points and polygons with Mapnik

texts:
I've tried to render a few maps in Mapnik and I have a little problem. I want to render a map that contains two layers: point and polygon. I've used this script:

import mapnik
rule = mapnik.Rule()
rule.symbols.append(mapnik.PolygonSymbolizer(mapnik.Color("darkgreen")))
rule.symbols.append(mapnik.PointSymbolizer())
style = mapnik.Style()
style.rules.append(rule)
layer = mapnik.Layer("mapLayer")
layer.datasource = mapnik.Shapefile(file="polska.shp")
layer.styles.append("mapStyle")
map = mapnik.Map(640, 480)
map.background = mapnik.Color("white")
map.append_style("mapStyle", style)
map.layers.append(layer)
map.zoom_all()
mapnik.render_to_file(map, "map.png", "png")

There is only one shapefile to render, but that is my problem. When I'm adding the 4th line rule.symbols.append(mapnik.PointSymbolizer()), I always get rendered polygons with additional points (centroids of each object I guess):

Why are there black points? Is there any possibility to render two layers (polygons and points) at the same time, without these centroids? What am I doing wrong?


texts:
I am trying to split a vector file using an attribute value and wants to name the each created output as the name in the attribute field. I am being able to split the vector layer but it is giving names by default to each of the file, i wants it to be specifically with the name in the attribute table. Suppose i have a field in attribute table title as NAMES and values in it as:
ABC
DEF
GHI
JKL
MNO
so i want it each to be split out and saved as NAMES_ABC, NAMES_DEF, NAMES_GHI, NAMES_JKL and so on
instead i am getting the names as NAMES_0, NAMES_1,NAMES_2, NAMES_3 and so on


texts: Default values from different layer in QGIS 3 - kind of spatial join

texts:
Is there a way to set the default value in QGIS 3 Properties - Attribute Forms based on the values of a different layer.
E.g. I set a point and the fields for e.g. country and height are filled out based on the position of the point and the values other layers.


texts:
I am trying to import information from a Bathymetric Attributed Grid (BAG) file into PostGIS using Python and psycopg2. I'm working, for now, on a Windows box with pre-compiled packages, so I don't have GDAL, osgeo, or ogr support for BAGs directly.

I can open and parse the BAG with h5py and ElementTree, no problem. I can add simple geo objects (e.g., bounding box as box2d) with no problems. I can't figure out how to add the BAG values. I'm trying to add the elevation and uncertainty arrays as two bands in a raster.

relevant, non-working code extract:

 elevation = hierarchy["/BAG_root/elevation"]
 width = len(elevation)
 height = len(elevation[0])
 uncertainty = hierarchy["/BAG_root/uncertainty"]
 el_array = elevation[0:][0:]
 un_array = uncertainty[0:][0:]

 update_string = "UPDATE %s SET el_un = ST_SetValues(ST_AddBand(ST_MakeEmptyRaster( %%s, %%s, %%s, %%s, %%s), ARRAY[ %%s, %%s ]::addbandarg[])) WHERE filename = %%s;" %(tablename, )
 print cursor.mogrify(update_string, (width, height, left, top, 1.00, el_array, un_array, basefn))

This gives a can't adapt type 'numpy.ndarray' error. I've tried several variations on this (without the slices, with tolist(), explicitly using numpy, etc., and all have had similar errors. I hope I'm missing something simple, due to my lack of PostGIS familiarity. For that matter, I may be able to copy the raster2pgsql output for another file type, suitably modified, but don't know how to encode the VALUES( '...'::raster) data value - is there some documentation for that I haven't been able to find?

I hope to soon move to a Linux box where I can add formats and recompile GDAL and other packages if needed, so answers using GDAL in Python or executables that are made with additional supported formats are welcome. I'd use raster2pgsql.exe, but it doesn't support the BAG format, at least as is. Is there an easy way to write out one of the formats it does accept by default?


texts: MULTIPOINT Text to list of POINT objects

texts:
I need to convert a MULTIPOINT string to a list of POINT objects. For example, suppose I have a query with ST_AsText() producing an output that is captured in my_points:

>>>print(my_points)
result={str}'MULTIPOINT(-1.0 1.2, 3.4 5.6, 7.8 9.0)'

I need result to be converted to a list of POINT objects so that I can use STRtree() to query if the points intersect a polygon.

I have tried the following:

[x.split(' ') for x in result[11:-1].split(',')]

but it only returns a list of strings:

(['-1.0', '1.2'], ['3.4', '4.5']...)

I need each element in each list to be floats so that I can convert them to Points.

What is the most efficient way to convert the MULTIPOINT string to an actual list of Points?

I am a novice to GIS and Python, not to programming in general.


texts: Developing workaround to arcpy.asShape truncating all values to 3 decimal points?

texts:
I've read a few questions\answers about converting a GeoJSON Shape into an ESRI Object.

Unfortunately, most of the answers pointed to decoding the dictionary's manually.

The problem:

The bug is that arcpy.asShape truncates all values to 3 decimal points and has been described at this site in arcpy.geometry __geo_interface__ and AsShape() function: loss of precision and holes. Naturally for most this is undesirable.

This is the solution i came up with:

from shapely.geometry import asShape
import arcpy

sr = arcpy.SpatialReference()
sr.factoryCode = 4326
sr.create()
arcpy.env.outputCoordinateSystem = sr

JSON = {'coordinates': [[[0.6049790000000712, 17.048294000000055],
                  [0.6068335000000502, 17.049157000000037],
                  [0.6067875000000527, 17.049668000000054]]],
 'type': 'Polygon'}

jsonobj = asShape(JSON)
arcshape = arcpy.FromWKT(jsonobj.to_wkt(),sr)

print (arcshape.WKT)
>>>u'MULTIPOLYGON (((0.60497900000007121 17.048294000000055, 0.60683350000005021 17.049157000000037, 0.60678750000005266 17.049668000000054, 0.60497900000007121 17.048294000000055)))'

Environment:

ArcGIS 10.1 SP1 for Desktop Quality Improvement Patch
ArcGIS 10.1 SP1 for Desktop Background Geoprocessing (64-bit) Quality Improvement Patch


texts: Filter the lowest / highest featureID in Geoserver CSS style rule

texts:
I am using GeoServer's CSS extension to style a Points layer which I publish as WMS. In my style, I would like to specify a rule that filters for the first and last features in my data set. Something like min(featureid) and max(featureid).

The documentation suggests how to filter for specific featureIDs or ranges of IDs. Consequently, I can filter for the first feature in my data set (called Points) and style it with

#Points.1 {
  stroke: #d73027;
  stroke-width: 6;
  stroke-linecap: round;
}

because the first feature will always have featureid=1. However, how can I determine the last featureID in my data set without hard-coding it since it might change over time?


texts: Arcgis 10: Point layer, polygon line layer, draw points to the nearest polygon point or other nearest line?

texts:
I have been using arcGIS for a while now but I often get frustrated with it :(

I have version 10.0 arcGIS and I'm have two layers

1) Point data (building locations)

2) a polygon shape file (outline of a continent)

A) I'm trying to find the nearest part of the polygon to the point data (nearest coastal point for the building) and draw a line between them forming a hypothetical pipeline.....

-> I have used spatial join and have a table of the nearest points but how to draw lines between those and the original points?

B) if these lines are close to each other or cross over, I want them to join up to form one line.

Any advice would be great here. All the options I tried are now either out of date or don't work anymore.


texts: Error using grass script v.lidar.edgedetection

texts:
I'm in trouble with the grass script v.lidar.edgedetection. When I try to filter an xyz format dsm the script starts but at the end of it an error message appears (see attached image).

Any idea where is the problem?

I tried to perform the same operation with a different file format (.las) but I have the same error.


texts: Mercator - create uniformly distributed dots over map

texts:
I have a map in Mercator projection. I know coordinates of corners. I want to render uniformly distributed dots to map.

for (float latDeg = startLatDeg; latDeg <= endLatDeg; latDeg += stepLatDeg){
        for (float lonDeg = startLonDeg; lonDeg <= endLonDeg; lonDeg += stepLonDeg){
            Pixel px = mercator->GetPixels(lonDeg, latDeg);
            map[px.x][px.y] = dot;
        }
    }

where start and end values are corners of map.

This code works for longitude, obviously, since it is uniform and non-distorted. How can I achieve uniformity in latitude? I know, that I cannot use fixed step, but somekind recalculate it based on last latitude. But I dont know how.


texts: ArcSDE: connecting via LDAP Authentication

texts:
I was wondering if anyone has a sample connection string for SQLPlus using LDAP Authentication. I am also trying to specify in the connection string the database I wish to use.

Thanks,
CDB


texts: Georeferencing raster with unknown highly non-rectilinear CRS in QGIS

texts: In QGIS, I'm trying to georeference an image where lines of latitude and longitude are curved / non-parallel (see below) and I have no information about the projection used. Can anyone recommend the transformation settings I should use and the CRS?


texts:
Sometimes the hardest part of an analysis is knowing what something is called.

What are some R packages, but more importantly, terminology I should be searching for to define contiguous city neighborhoods based on home price changes - basically a topographic map treating these price changes as elevations (and ignoring how real estate people have defined neighborhoods).

I'm not sure if this is called "natural neighborhood delination" or "spatially restricted clustering" or just "grouping analysis".


texts: QGIS - Create equal length polygon along polyline

texts:
How can I create a series of polygons along a line? I need to assess the density of a number of targets along a line within a corridor.

I have tried the following:

creating a buffer,
using the plugin station lines to create perpendicular lines,
Using the Digitizing Tools Plugin to split the buffer with the
station lines

This only splits the polygon one line at a time, the lines is 20+km long and I would like 25m sections making manual repetitions impractical.The line also has complex curves.

Is there a way to create such a grid?


texts: How to Delete Fields in bulk using Modelbuilder?

texts: Is there a way to simultaneously delete fields in Modelbuilder? I want to be able to remove excess fields from an iterative join of near tables. All of the fields I want to delete and keep are similarly named (ex: NEAR_FID). I have 40+ copies of that variables named this way in one table.


texts: How to get the affine matrix between 2 projection using pyproj?

texts:
I'm struggling getting the affine matrix (or crs_transform) using pyproj. The Transformerobject is full of functions to actually transform coordinates but I just want to check the matrix:
I'm stuck here:
import pyproj

toto = pyproj.transformer.Transformer.from_crs("EPSG:3857", "EPSG:4326")
toto.to_json()


texts:
I have a map with some markers. On clicking a marker, a popup would appear and on clicking the map, any visible popup would be hidden. I want to hide any currently visible popup upon clicking a marker as well, and introduced a line of code in the following event listeners:

marker.events.register("click", marker, function(event) {
 map.events.triggerEvent("click", map); // this line of code added
 popup.show();
});
map.events.register("click", map, function(event) {
 popup.hide();
});

However, I get an error upon clicking of a marker with the error given as "px is undefined" on line 1981 in Openlayers.debug.js file, which is the Openlayers.Pixel class with distanceTo method. I am not sure what is wrong and can't seem to make sense of it.

Can anyone help me? Thanks.


texts: Is there a way to re-order how layers render in Mapnik?

texts: I'm creating Mapnik stylesheets and I am trying to edit the order layers are rendered. I tried just moving the XML code, but that didn't seem to work. Are there specific steps to make some features render above others?


texts: 3d interpolation between two xyz coordinates

texts:
I need to find the x and y coordinate on a known z coordinate based on two known xyz coordinates.  In more detail: I have a csv-file with x y and z values and I need to find the place where the 0, z value is crossed.

The data looks like this:

z,x,y
0.822,174269.8332,608470.8717
0.692,174269.5191,608480.8668
0.502,174269.205,608490.8618
0.382,174268.8909,608500.8569
0.272,174268.5768,608510.852
0.172,174268.2627,608520.847
0.022,174267.9486,608530.8421
-0.262,174267.6345,608540.8372
-0.772,174267.3203,608550.8322
-1.052,174267.0062,608560.8273
-0.792,174266.6921,608570.8224
-0.382,174266.378,608580.8174

The (very basic...) illustration below illustrates what I mean. The green and the red lines are the known points and the blue line between those points is a slope that I could calculate I suspect. The yellow point with the circle around t is the x and y coordinate that I have to find:

I have been looking at scipy.interpolate.interp2d but that returns the error:

Traceback (most recent call last):
  File "D:\ARNO\JarkusPython\ARNO\Create_Data\createStraightData.py", line 50, in <module>
    some = scipy.interpolate.interp2d(xx, yy, zz)
  File "C:\Python27\ArcGIS10.5\lib\site-packages\scipy\interpolate\interpolate.py", line 229, in __init__
    self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
  File "C:\Python27\ArcGIS10.5\lib\site-packages\scipy\interpolate\fitpack.py", line 929, in bisplrep
    raise TypeError('m >= (kx+1)(ky+1) must hold')
TypeError: m >= (kx+1)(ky+1) must hold

Which means (I believe) something as that there's not enough coordinates. Also, I think I have to feed the function the x and y and not the z value. Is this the mathematical problem that I run into?

And does anybody know a solution to this?


texts: SAGA error (Segmentation fault) with running QGIS

texts:
--Running with configuration is Mac OS X 10.11.6 (El Capitan), QGIS 2.14.7-Essen, SAGA 2.2.3--
I'm hoping that someone knows what the fix is for the failure I'm encountering in the SAGA tool.

I'm following the instructions in the QGIS Training Manual (http://docs.qgis.org/2.14/en/docs/training_manual/vector_analysis/spatial_statistics.html), Section 7.4.11 for running SAGA's Minimum Distance Analysis tool. Instead of a result, the GUI reports a failure in the SAGA run with a Segmentation Fault 11.  The GUI states that the Minimum Distance Analysis layer was not correctly generated. Just to note in case the projection is pertinent, the random points used in the analysis are in the WGS 84/UTM zone 34S projection.

/Users/xxxxx/.qgis2//processing/saga_batch_job.sh: line 3: 1625 Segmentation fault: 11 saga_cmd statistics_points "Minimum Distance Analysis" -POINTS "/var/folders/w5/0qqndryj29587bhbbt_wxx1w0000gn/T/processing8e4aa8219d5c48cfbc84b90c68593311/1485958248.417.shp" -TABLE "/Users/xxxxx/Downloads/junk.csv"
Converting outputs
Loading resulting layers

The following layers were not correctly generated.
Minimum Distance Analysis
You can check the log messages to find more information about the execution of the algorithm

I reran the tool via the command line to see if I could find more information.  (See below.)  The only additional clue I see is something about "Load shapes".

/Applications/QGIS.app/Contents/MacOS/bin/saga_cmd statistics_points "Minimum Distance Analysis" -POINTS "/var/folders/w5/0qqndryj29587bhbbt_wxx1w0000gn/T/processing8e4aa8219d5c48cfbc84b90c68593311/1485958339.968.shp" -TABLE "/Users/xxxxx/Downloads/junk.csv"
____________________________

   #####   ##   #####    ##
  ###     ###  ##       ###
   ###   # ## ##  #### # ##
    ### ##### ##    # #####
 ##### #   ##  ##### #   ##
____________________________

SAGA Version: 2.2.3

____________________________
library path:   /Applications/QGIS.app/Contents/MacOS/lib/saga/
library name:   libstatistics_points
library     :   Points
tool        :   Minimum Distance Analysis
author      :   O.Conrad (c) 2010
____________________________

Load shapes: /var/folders/w5/0qqndryj29587bhbbt_wxx1w0000gn/T/processing8e4aa8219d5c48cfbc84b90c68593311/1485958339.968.shp...

Segmentation fault: 11


texts: Interpolation from isolines vs Thiessen (Voronoi) polygons interpolation

texts:
I'm currently studying GIS and have come across a question in an assignment that I just can't answer. I've been googling my brains out for the last few hours and got nowhere... here goes...

Identify two (2) similarities and two (2) differences between (i) interpolation from isolines and (ii) Thiessen (Voronoi) polygons interpolation. Explain your answer.

Any ideas or know somewhere to point me where I might find the info?


texts:
I use MapInfo and have a lot of data in .tab format that I would like to put into a database in order to access it through a website. I would like to continue using MapInfo as my GIS programme and export to the web database with the minimum of fuss.

What would be a good type of database for this?

I was considering MySQL as I have used it a lot before and it has native spatial data types. However I have read that MapInfo can't connect directly to MySQL and that its spatial data support is a bit patchy.

Alot of discussion on the subject is several years old though, so can anyone confirm or deny this?

If not MySQL, what is a good choice? PostgreSQL/PostGIS?


texts: can grass gis module be used in parallel mode

texts: i want to use my customized grass gis(6.4.3) module for separate bands of a sensor in a parallel mode. is it possible? i have created separate location and if one band is executed ,for the second band message appears as " not found".I searched in the grass manual and seen "GIS_LOCK" variable. can any one suggest me how to use this for my executable in parallel mode.


texts: Converting miles above ground to zoom level?

texts:
I've got a spec that details what features are to appear on a map at certain number of miles above the earth.  So for example "Render this point when closer than 400 miles above the earth".

How do I convert that to a scale? I'm using a standard zxy base map so I need to change the spec to "render this point when zoom level greater than 6"

Is there a standard calculation for this?


texts: C-based implementation of Djikstra's or Fast Marching for QGIS/Python?

texts:
I need an algorithm to find least cost paths between points on a discrete grid,

allowing for the direction of transport to affect the cost of transport,
allowing for different modes of transportation which connect
in some directions and not in others, and
allowing for "transshipment" costs for switching modes of transport.

For this particular application, I need the algorithm to be blazing fast as it    will be executing many, many times as part of a larger optimization routine.

For MATLAB, I have found a C-based (MEX) implementation of the Fast Marching Method, which with a bit of modification I have succeeded in generalizing to a Djiktra's Algorithm implementation with the features I need. This is fast enough, but it would be cleaner and more convenient to be able to run my entire optimization routine out of QGIS, as this is where a lot of my data needs to be processed anyway.

What are my options? What do you recommend?


texts: Connecting isolines (i.e.points of equal value) in 3D scene using ArcGIS for Desktop?

texts:
I have a point shp theme with fields Lat, Long, Altitude (m), and Number.
The Number field – is a numerical parameter of each point (Lat, Long).

The Altitude field – is the parameter of elevation of each point (Lat, Long) above the surface.

I made a 3d map in ArcScene of ArcGIS10 where I represented each point as a symbol of different size (according to the “Number” field), and then – elevated each point above the surface (according to the “Altitude” field).

Is it possible to connect by isolines (or something like that) all  points of the equal value of the “Number” field (i.e. all points of the same size) in this 3d map?


texts: Using local tiles created in Geowebcache with OpenLayers

texts:
I am trying to setup an Openlayers page that uses a local tile cache that I have created using Geowebcache. But the HTML page and the tiles will be stored in an Amazon S3 bucket
So far I have setup GeoServer, created a custom gridset and seeded all the zoom levels.
I am having problems writing an OpenLayers HTML page that reads those local tiles.
My tile cache structure looks like this
zoom00/0_0/0_0.png
zoom00/0_0/0_1.png
etc

In my OpenLayers I have the following
var layerExtent = new OpenLayers.Bounds( 0,  0, 700000 , 1344000);
map = new OpenLayers.Map( 'map', {'restrictedExtent': layerExtent} );
layer = new OpenLayers.Layer.XYZ( "Data Layer",
"https://url_goes_here/${z}/${y}/${x}", {sphericalMercator: false} );
map.addLayer(layer);
map.zoomToExtent(map.restrictedExtent);

After some google I found this link to a Dev example that seems to have a geowebcache as an openlayers.layer
http://dev.openlayers.org/sandbox/tschaub/gwc/examples/geowebcache.html
But when I look at the API docs there isnt such a layer? So maybe that it doesn't exist.
I then found an OSM Local Example here
http://wiki.openstreetmap.org/wiki/OpenLayers_Local_Tiles_Example
Which uses
var newLayer = new OpenLayers.Layer.OSM("Local Tiles", "tiles/${z}/${x}/${y}.png", {numZoomLevels: 19, alpha: true, isBaseLayer: false});
But again their PNG files see to have a different naming convention than the ones created by Geowebcache.
Is there a way to change the Geowebcache naming convention or is there something I can do in the OpenLayers code?
thanks for any help
UPDATE
Thanks to CHenderson I have a working solution.
using his option of just downloading the GeoWebCache.js and including it in my HTML I created the following:
<script src="http://openlayers.org/api/OpenLayers.js"></script>
<script src="geowebcache.js" type="text/javascript"></script>
<script type="text/javascript">
     var map;
        function init() {

            var world = new OpenLayers.Bounds(
                0, 0, 700000, 1344000
            );

            map = new OpenLayers.Map({
                div: "map",
                projection: new OpenLayers.Projection("osgb_zoom"),
                numZoomLevels: 11,
                resolutions: [896,448,224,112,56,28,14,7,4,2,1],
                maxExtent: world,
                layers: [
                    new OpenLayers.Layer.GeoWebCache({
                        url: [
                            "https://url_to_tile_cache/"
                        ],

                        buffer: 1
                    })
                ],
                center: new OpenLayers.LonLat(400000, 500000),
                zoom: 0
            });
        }

        OpenLayers.Util.onImageLoadError = function() {
            this.src = "error.png";
            this.style.display = "";
        };

</script>

Even requesting from just one tile cache is impressively quick.


texts: How can I bring rasters into R using a previously defined year in the name

texts:
This seems like A simple question but I'm not sure how to phrase it or ask.
I want to bring in to r many rasters from many years and instead of typing them all out, I want to be able to define the year and bring in rasters by incorporating that defined year. This way I can copy and paste my code and instead of changing the year in every raster name, I can just change it once where I define the year.
Here is an example. I defined year as 2012 and in r.12 I tried using it to bring in the raster but this didn't work. r.1 is an example of how the rasters are all named.
year <- "2012"

r.12 <- raster("pr_total_mm_CRU_TS40_historical_12_"+year+".tif")
r.1 <- raster("pr_total_mm_CRU_TS40_historical_01_2012.tif")


texts: How to get latitude and longitude for a matrix data using specific projection in R?

texts:
I have matrix of temperature data and I would like to get lat and long for this data set in a particular projection

e.g.

PROJCS["DHDN_3_Degree_Gauss_Zone_3",GEOGCS["GCS_Deutsches_Hauptdreiecksnetz",DATUM["D_Deutsches_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128]]

Does anyone know how do this in R?

The matrix data is like:

NCOLS 654   654
NROWS 866   866
XLLCORNER   3280414.    711633467
YLLCORNER 5237500.  62890625
CELLSIZE 1000
NODATA_VALUE -999
-999    -999    -999    -999    -999    -999    -999    -999    -999    -999 -999 15 18 12 13 20 15 16 18 21


texts: Relationship of QGIS to its sister programs?

texts:
I'v e been using QGIS 2.4 for a number of months now and we are using it a bit at work.

I noted that when I did the initial install a number of new icons appeared on my desktop namely
MSys
QGIS Desktop 2.4.0
QGIS Browser
Grass
OSGEOWShell
Saga GIS (2.0.8)

If I need to roll this out to other machines I would prefer just to have the minimum footprint. I can see that I probably need the browser as that looks like it is within QGIS Desktop but Grass / OSGeowShell and Saga as well.

I'd be interested if anyone could indicate any interdependencies between the programs and what is the minimum install. I suspect its Desktop and Browser?


texts: I am aware how to rasterize a 2D line, but if I have two coordinate triplicates (x1,y1,z1) and (x2,y2,z2), how can I efficiently rasterize the resulting connecting line such that the raster cell values represent the proper slope of the original points? I am willing to explore any solutions through ArcGIS, Arc Workstation, GDAL/OGR, SAGA GIS, etc.


texts: Is the Intersection Operation Commutative?

texts:
In map overlay, is intersection of map layers a commutative operation ie layer_A intersection layer_B = layer_B intersection layer_A?
If layer A is a square and layer B is a triangle. In this case, (A intersectn B) == (B intersectn A) ?


texts: How to set the number of vertices by shape digitizing tool in QGIS 3?

texts:
Can we set the number of vertices of the created circle or ellipse.
(e.g. a circle/ellipse with 24 or 36 vertices)
I can use "add regular polygon and set the point #" for circle, but i don't know how to set the ellipse.


texts: Label expression in ArcGIS gives TypeError coercing to Unicode need string or buffer NoneType found?

texts:
I try to build stacked label expression in ArcGis  10.3.1 with this expression:

[name3] + '\n'+ [results]

name3 and results are two columns in the attribute table. But I get this error:

typeerror coercing to unicode need string or buffer nonetype found

I read Building label expressions - Expression examples but didn't find any solution.

This is the attribute table:


texts: Can I import raster files from GRASS GIS in QGIS versions later than 1.8.0

texts: I have a process that I have used for years where I create .ascii files with a Python process, import them to GRASS GIS 6.4.3 and do some raster calculations (r.mapcalc).  This works great but then for mapping in QGIS, I have to go back to version 1.8.0 (I'm glad I have saved this old version!) to create output maps of my calculated raster maps.  Am I missing something or has QGIS left compatibility with GRASS rasters behind?


texts: Getting Annotation selection using ArcObjects?

texts:
I have annotations stored in mxd document, when I select some annotation, I tried to get the selection through graphics selection, so I always get the rectangle of the dataframe (because I'm working on a axPageLayoutControl), and when I use feature selection I get nothing; so  none of them return the annotation selection.

How could I get the annotation selection in ArcObjects?

I'm using ArcObjects 10.1.


texts: How do I project my WMS tilelayer in leaflet if I want it to display as an overlay on my basemap if they have different EPSG?

texts: So my basemap is a world Mercator projection and in wfs format and I'm trying to put a WMS overlay in a different projection over it? Leaflet's API didn't give me much info on changing the epsg so is it even possible?


texts: I am having a problem.  I follow the steps to try to georeferenced, but I see that the JPEG is much bigger than the shapefile when you add both to the table of contents.  I try changing the extent on one, but it doesn't make it match the other.  I tried setting extent to data frame, but it still does not change it.  I am trying to georeferenced a JPEG map of the U.S. territorial acquisition to a shapefile downloaded of the U.S.  The JPEG map shows up so much bigger than the shapefile that I cannot even see the shapefile to make any control points or anything.  Here is a snapshot of the problem.  It is just about impossible to see the U.S. shapefile in the background, but if you look carefully at the top left corner you can see it.  There is no way for me to georeference this image because the image is way too big compared to the shapefile, and nothing seems to change if I try to change the extent.


texts: Zonal statistics returns empty on Landsat 8 raster data in QGIS

texts:
I have derived surface temperature from the Landsat 8 thermal band 10.

My vector file is the following: Small Areas City of Melbourne

I want to calculate the mean, median and maximum temperature of the 14 different small areas that are within my spatial extend. I have already clipped the temperature raster file to the outline of the small areas vector layer.

Now I have tried Raster > Zonal statistics which simply does not return anything
and I have tried the Zonal Statistics from the plugins, which returns empty values.

Can someone help me out?


texts: ST_MAKEVALID: ERROR: Geometry type (MultiPolygon) does not match column type (Polygon)

texts:
I'm using Postgres 9.4 and PostGIS 2.1. I have a table that looks like this:

 wkb_geometry | geometry(Polygon,4326) |
 id           | integer                |
 county       | character varying      |

I've learned that the table contains some invalid polygons (e.g. when I try to do ST_UNION I get errors), so I'm trying to fix them:

UPDATE mytable SET wkb_geometry=ST_makevalid(wkb_geometry) WHERE id='1';

But this gives me an error:

ERROR:  Geometry type (MultiPolygon) does not match column type (Polygon)

What should I do about this?


texts: gdaldem hillshade stop working with new projection?

texts:
Parallel to gdalwarp `-te -5.8 41 10 51.5 -ts 1980 0` not working with mercator? , my hillshade works with lat/long input :

gdaldem hillshade input.tmp.tif shadedrelief.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges

But stopped to work when input is in mercator :

gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:3857  input.tmp.tif reproj.tmp.tif #mercator
gdaldem hillshade reproj.tmp.tif shadedrelief.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges

While lat/long is in degree, mercator is in meters, I believe both z scales to be in m which requires -s 111120 (requires -s 370400 if in feet!). Also, what do I do wrong ?


texts: Entering multiple image on single page using ArcGIS Desktop?

texts:
I have 10 dynamic images.

I want to print these images (10) on a single layout page.

How to do this using ArcGIS 9.2?


texts:
I need to prepare a 1m contour map for an area with LiDAR data (XYZ files, in a 1m grid). The files are a little large and there are 54 of them. I merged the files to a single XYZ file (which turned out to be a 12GB file BTW) and imported it as a vector layer using delimited-text option.
I was hoping to use TIN interpolation and then to use contour extraction. But, for some reason, TIN interpolation is stuck at 0% even after a couple of hours.
(with pixel size set to 1.0)
I apologize if this question might be silly or I'm violating some basic principle, but I'm totally new to this. I'm not even sure if import as vector --> TIN --> extract contours is the correct way to go about things, but that's the way I saw mentioned in most places.
I tried importing individual XYZ tile, convert to TIN, then extract contours first, and it worked. But, then each tile didn't match the adjacent one.
Also, I made the merged file with this Python program, which I hope is correct (it did managed to merge the files successfully):
import os
ABS_PATH = os.path.abspath(__file__)
BASE_DIR = os.path.dirname(ABS_PATH)
DATA_DIR = os.path.join(BASE_DIR, 'xyz_tiles')
files = os.listdir(DATA_DIR)
with open('merged.txt', 'w') as merged:
    for file in files:
        with open(os.path.join(DATA_DIR, file)) as f:
            for line in f:
                merged.write(line)


texts: Exporting maps at high resolution for a specified area

texts:
I am using ArcGIS 10.1 and am very new to it. I have Bing Maps aerial imagery of a particular road intersection and I want to export the map as a TIFF at a 1:1000 scale. When I zoom in and then use Export map feature only the area which I can see on screen gets exported. If I zoom out into the whole study area I get a very bad resolution TIFF. I know the top, bottom, left and right extent of my area in meters. Now, my question:

Is it possible to get multiple images (at 1:1000 scale) of the area - like tiles -  that can be joined to get the whole area map? A detailed procedure description highly appreciated.


texts: What's the difference between ArcGIS Engine and ArcGIS for Server?

texts:
Can anyone provide a succinct answer as to the difference between ArcGIS Engine and ArcGIS for Server?  I am working on a justification for Server and I want to make sure I don't miss any key points.

A couple of similar and related questions which did not prove useful to me are:

Difference between ArcSDE and ArcGIS Server?
What are the differences between the various ArcGIS Server options?


texts: Create junction raster from stream and road rasters

texts:
I have converted two vector polyline datasets to raster format. One represents streams, and one represents roads. I need to create a third raster layer that contains cells where there are stream crossings.

I was going to use a simple boolean AND function to create a raster of stream-road junctions. However, at about 10% of the junctions, the raster lines intersect, but do not overlap. See the figure below.

Is there a tool or established method that edits the road raster to make sure that all road crossings actually have overlap with raster streams dataset?


texts: When would a PostGIS geometry compare using ~= return true when ST_Equals returns false?

texts:
I have a couple of featureclasses in a 10.2 geodatabase/PostgreSQL 9.2.2 database and I've been struggling with the FeatureCompare tool--it is telling me geometries are different on some features where I know them to be at least spatially equal.

Luckily, I have an ID and I could join between the 2 tables in PGAdmin and compare geometry.  The ~= opreator and the ST_Equals both return true for every row.  Then I modified 1 feature by moving 1 vertex 800 ft to the South.

When I run my query, ST_Equals returns false which is what I expect.  But ~= still returns true??? Why is that?

select t1.shape,t2.shape ,(t1.shape ~= t2.shape) as exactequal , ST_Equals(t1.shape , t2.shape) spatialequal, (t1.shape = t2.shape) bndequal
from  pg.Test1 t1
inner join pt.Bnd2 t2
on t1.polygon_id=t2.polygon_id
where t1.objectid = 6389

I trying to figure out if this is a bug or if I'm misunderstanding the operators.

"PostgreSQL 9.2.2, compiled by Visual C++ build 1600, 64-bit"
"POSTGIS="2.0.4" GEOS="3.4.2-CAPI-1.8.2 r0" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24"


texts: How to calculate the intersection within the grid and saving to the grid

texts:
Here is my example:
http://i4.tietuku.com/dab1aaa0831739c8.png

(1) I creat the grid line layer using MMQGIS plugin
(2) I use line intersection function to get the intersection points.

I want to calculate the distance sum within each grid and save the grid line layer into raster data.

I'm using QGIS 2.8.3.

Add

Sorry for my bad expression.
The polyline represent the road network.
And I want to calculate the distance of road network within the grid which present the road length in this area.


texts: Geoserver/Geofence - Trigger Authentication

texts:
Geoserver Catalog MIXED mode in Data Security section is perfect for my purposes since I want to trigger an alert for authentication whenever a specific layer is requested.

Since I need to create Geoserver Users using a REST API from another server I need to use Geofence.

When defining rules in Geofence, there's an option to modify the catalog mode for a specific layer.

I expected this "authentication alert" to appear when the catalog mode in Geofence rules is defined as MIXED. But this is not happening, instead, I receive this error:

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://localhost:8080/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException code="LayerNotDefined" locator="layers">
      Could not find layer workspace:layer
</ServiceException></ServiceExceptionReport>

Interestingly, If I navigate to the dev tools and check for the request, I'm getting some kind of token in the body:

So my questions:

Is this the expected behaviour?
Or i'm doing something wrong?
What is this token I'm receiving?


texts: Is there a faster way to get unique values for a given FeatureClass in ArcObject?

texts:
I'm trying to count rows by unique values for a given FeatureClass. I'm using DataStatistics to get the unique values for a specific column for a given FeatureClass.

This seems to be (very) slow with big shapefiles (> 200k of features) - is there a faster way to determine unique values?

The second step of the process (using QueryFilter in combination with FeatureClass.FeatureCount()) seems fast enough although I would be glad if im pointed to a better solution.

public static Dictionary<string, int> CountUniquesStatistic(ITable table, string fldName)
{
    Dictionary<string, int> uniqueValuesDictionary = new Dictionary<string, int>(StringComparer.CurrentCultureIgnoreCase);

    DataStatistics dataStatistics = new DataStatistics();
    dataStatistics.Field = fldName;
    dataStatistics.Cursor = table.Search(null, false);

    IEnumerator uniqueValues = dataStatistics.UniqueValues;

    while (uniqueValues.MoveNext())
    {
        object current = uniqueValues.Current;
        uniqueValuesDictionary[current.ToString()] = 0;
    }

    foreach (string key in uniqueValuesDictionary.Keys.ToList())
    {
        IQueryFilter queryFilter = new QueryFilter();
        queryFilter.WhereClause = string.Format("{0} = '{1}'", fldName, key);

        uniqueValuesDictionary[key] = table.RowCount(queryFilter);
    }

    return uniqueValuesDictionary;
}


texts:
I want to set up a layer selector with CartoDB like what is found in Usage Examples here. I viewed the source code and changed it for the table I want to use. When I open the html file and click on a different layer, I get a 404 error. The URL for the error states:

({"errors":["Error: could not fetch source tables: column \"undefined\" does not exist"]});

But I can see and query my table in CartoDB just fine. Below is my createSelector function. I'm not trying to do anything more than what is in the example I linked above. Does anyone have an idea on what I am missing? To be clear, I'm not sure the error above is related to a problem with my code, but it could be! I have specified my CartoDB user name and my json file.

Collisions

Pedestrian 2013

Pedestrian 2012

BIKE 2013

Bike 2012


texts:
Is there anyway to get the Identify Features button in QGIS to automatically select a layer and provide the information without having to highlight it first in the layers TOC?

MapInfo and ArcGIS automatically highlight, but can't see that QGIS does.


texts: Elevation info for multiple lat long coordinates

texts:
I have a list of 100 lat long pairings that I wish to get the elevation data for.

I am familiar with Google Earth and I can plot each of these points but when I check the properties there is no elevation info other than the info at the bottom tool bar which would require me to manually write the elevation down.

Is there another way that I can query all 100 coordinates to get the elevation data?


texts: Problem compiling ui to python for a QGIS 3 plugin

texts:
I've created a .ui using QT Designer for a QGIS 3.10 plugin called 'bhlayers' which has a docking widget to the right. When I used the: pb_tool deploy command it didn't seem to translate the .ui to the bhlayers.py file as I expected. It hasn't set up the Qt buttons, checkboxes, horizontal layouts, labels, lineEdits etc.

Is it not supposed to?  If it isn't how do you relate the .ui objects with the python code that detects interactions with the interface. Example:

if dbBox.isChecked() == True:  ..

So I then tried instead to use pyuic5:

pyuic5.bat -x bhlayers_dockwidget_base.ui -o bhlayers.py

This time the bhlayers.py has all the Qt objects created, buttons, checkboxes, etc.

However it looks like pyuic5 is compiled in a completely different format to pb_tool and is not compatible (at least without a lot of editing) with how QGIS plugins work.

I've created QT applications before but manually without the QT Designer. I wonder whether I should have stuck to that - as I am a bit lost now.

Does anyone know what I'm doing wrong? Or what wrong assumptions I'm making?


texts: Is there an open source equivalent to ArcGIS "collect events"?

texts:
I am looking for an equivalent tool to ArcGIS "collect events" (within the Spatial Statistics -> Utilities toolbox).  I wish to aggregate coincident point data in a space where many points lie in exactly the same place making visualisation difficult.

For example, to go from:

to:

thus providing a visual guide on density on points before proceeding to other forms of spatial analysis.

Is it possible to perform this on a different platform, preferably QGIS (or the associated PostGIS spatial db)?


texts: Auto populate field based on another field in Esri Collector

texts: I'm looking to have an automatic conversion worked into an extra couple fields that will convert the latitude and longitude taken when collecting a point. Currently trying to convert the latitude and longitude to state plane coordinates. I've downloaded a toolbox that creates the latitude and longitude automatically from my Trimble gps receiver.


texts: Permission denied when trying to Initialize google earth engine

texts:
I was following the google earth engine documentation here. I added all the files in the root of my project as per the documentation. These files are the config.py and the privatekey.json. I then wrote the following code:
import ee
service_account = '<service-acc>@<project-id>.iam.gserviceaccount.com'
credentials = ee.ServiceAccountCredentials(service_account, 'privatekey.json')
ee.Initialize(credentials)

When I run this code, I got this error:
Traceback (most recent call last):
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/ee/data.py", line 345, in _execute_cloud_call
    return call.execute(num_retries=num_retries)
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/googleapiclient/_helpers.py", line 134, in positional_wrapper
    return wrapped(*args, **kwargs)
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/googleapiclient/http.py", line 907, in execute
    raise HttpError(resp, content, uri=self.uri)
googleapiclient.errors.HttpError: <HttpError 403 when requesting https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/algorithms?prettyPrint=false&alt=json returned "Permission denied.">

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/ee/__init__.py", line 125, in Initialize
    ApiFunction.initialize()
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/ee/apifunction.py", line 154, in initialize
    signatures = data.getAlgorithms()
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/ee/data.py", line 1052, in getAlgorithms
    return _cloud_api_utils.convert_algorithms(_execute_cloud_call(call))
  File "/home/thuha/miniconda3/envs/prj/lib/python3.8/site-packages/ee/data.py", line 347, in _execute_cloud_call
    raise _translate_cloud_exception(e)
ee.ee_exception.EEException: Permission denied.

I need to know what makes that error appear when I want to initialize.
Google documentation


texts: Difference between Image Analysis Window and other associated tools?

texts:
I've been unable to find some concrete ESRI documentation on this.

What are the differences between the Image Analysis Window tools and the corresponding tools in other toolboxes? How come the Image Analysis is so much faster?


texts: More efficient way to pan and mark locations viewed in ArcMap10

texts: I'm reviewing a large dataset and I would like to be able to pan through it with the keyboard arrows, and somehow mark the areas that I've viewed. Is there a way do this more efficiently than dragging the mouse across the screen for every pan?


texts: Apply mann-kendall analysis on stacked raster in R as a function

texts:
I want to run a mann-kendall trend analysis with pre-whitening on a stacked raster and then save the resulting tau, sen-slope and p_value.  I am currently doing this by looping through every pixel but it is incredibly slow so I am wondering if there is a better approach to this? I am using the modifiedmk package for the trend assessment. I would be ok with any package that applies pre-whitening, and returns a sen slope, a tau value and a p-value though.

This is my current code:

library(Kendall)
library(zyp)
library(gimms)
library(raster)
library(rgdal)
library(modifiedmk)

#create a fake raster
r = raster(nrow = 50, ncol = 50)
r[] = rnorm(n = ncell(r))

#make raster stack
stacked = stack(x = c(r, r*2, r))
stacked2 = stack(x = c(r, r*2, r))
stacked3 = stack(x = c(r, r*2, r))

#get a depth of 9
final_stacked = stack(stacked, stacked2, stacked3)

#some empty data to repopulate
coef = raster(nrow = 50, ncol = 50)
p = raster(nrow = 50, ncol = 50)
t = raster(nrow = 50, ncol = 50)

#loop through every pixel stack and runb a mann kendall test
for (r in 1:nrow(coef)){
  for (c in 1:ncol(coef)){

    #one pixel stack
    pixel = data.frame(final_stacked[r, c][1,])
    names(pixel) = 'Values'

    #run mann kendall
    result = pwmk(pixel$Values)

    #save the tau, p-value and sen slope
    tau = result[[7]]
    p_value = result[[4]]
    sens = result[[3]]

    #save the value to new rasters
    coef[r,c] = sens
    p[r,c] = p_value
    t[r,c] = tau
  }

}

I also tried this which should be faster than looping if it worked:

tsfun <- function(x) {
  if(all(is.na(x))){
    return(NaN)
  } else {
    r <- data.frame(pwmk(x))
    names(r) = c("Z", "Sen", "Old_Sen", "p", "s", "var(s)", "Tau")
    a <- r$Sen
    return(a)
  }
}

raster.mkt <- calc(final_stacked, fun=tsfun, na.rm = TRUE)

but it returns:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function. Perhaps add '...' or 'na.rm' to the function arguments?


texts: How to add longitude / latitude to a world map in MapInfo professional?

texts:
I am a novice in MapInfo Professional. I have imported a shapefile of the world map into MapInfo. It's projected as Lambert azimuthal equal-area projection (ETRS89). Is it possible to show the longitude / latitude as well?

After what I have read in the forums, I am guessing I could do something with Grid maker or grid tools?

I can't find any specified steps for this requirement.


texts: How to import/export geographic information (gml or csv file) to geodatabase?

texts:
I am searching for C#/Java libraries that can import and export geographic information (gml/csv file formats) from and to geodatabases.

If not yet available, is there other way to do this task?

Any guidance on the right direction is appreciated.


texts: Geoserver WMS Cascading to MapXtreme, PNG not transparent

texts:
Geoserver allows to proxy WMS GetMap and GetFeatureInfo Requests to other WMS-Servers. Overall this works reliably and great.

However, I'm having trouble with a new service I'm integrating. No matter what I do, resulting PNG's always have a white background instead of transparent. The exact same Request pointed directly at the MapXtreme WMS Service works as it should, meaning the resulting PNG is transparent.

Everything is configured as it should, as far as I know. Moreover, it is configured exactly the same as an other working WMS (I don't know which technology is used there).

The thing that stood out for me in the log file was the BgColor parmeter:

2018-03-22 09:22:51,878 INFO [geoserver.wms] -
Request: getServiceInfo
2018-03-22 09:22:51,880 DEBUG [geoserver.ows] - Getting layers and styles from LAYERS and STYLES
2018-03-22 09:22:51,880 DEBUG [geoserver.ows] - establishing null style for myws:mylayer
2018-03-22 09:22:51,881 DEBUG [geoserver.wms] - setting up map
2018-03-22 09:22:51,881 DEBUG [wms.map] - setting up 256x256 image
2018-03-22 09:22:52,326 INFO [geoserver.wms] -
Request: getMap
    Angle = 0.0
    BaseUrl = http://my.wms.server:80/
    Bbox = SRSEnvelope[2717396.2137828125 : 2722683.5524514774, 1248403.6496270949 : 1255733.3898753526]
    BgColor = java.awt.Color[r=255,g=255,b=255]

What did I miss, why am I not getting transparent PNG's?

Edit:

WMS Parameters are the following, sent to both services. One is transparent, one isn't:

SERVICE: WMS
VERSION: 1.1.1
REQUEST: GetMap
FORMAT: image/png
TRANSPARENT: true
STYLES:
LAYERS: my_layer
SRS: EPSG:2056
WIDTH: 256
HEIGHT: 256
BBOX: 2718121,1248957,2723408,1256286

URL encoded as: SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=my_layer&SRS=EPSG%3A2056&WIDTH=256&HEIGHT=256&BBOX=2718121%2C1248957%2C2723408%2C1256286

Edit:

My colleague found out, MapXtreme produces RGB PNGs with NoData values instead of RGBA. Geoserver can not work with these, thus producing a PNG with a white background. There are some issues in GeoWebCache related to this. We are currently working to get this fixed in Geoserver.
The other, easier option would be to configure MapXtreme WMS to produce the right output format for PNG. If this is at all possible, and if so, how?


texts: Getting wrong value for area of geometry using Spatialite?

texts:
I am trying to implement a measuring tool on apple maps for IPAD. As a part of that i collect all the touch points on the map into an array and constructing a geometry object from that and then calculating the length and area of the geometry using Spatialite.

Now my problem is that i am getting the length calculated correctly but the area is not calculated correctly

Below is the query for calculating length of geometry

SELECT Length(Transform(GeomFromText('POLYGON((32.653995 39.754341,32.654217 39.754127,32.653842 39.753991,32.653647 39.754207,32.653995 39.754341))',4326),3857))

Expected value is around 131 and Query result is 132.

Below is the query for calculating area of geometry.

SELECT Area(Transform(GeomFromText('POLYGON((32.653995 39.754341,32.654217 39.754127,32.653842 39.753991,32.653647 39.754207,32.653995 39.754341))',4326),3857))

Expected value is around 1020 Sqm but Query result is 1706 Sqm.

Note: all points are geo referenced to Turkey region

Why is this happening?


texts: Extracting roads out of Landsat raster images

texts: For a study into the expansion of road networks in the rainforest I am trying to extract roads out of Landsat images. We already have sharp and cloud-free composites on which the roads are clearly visible by eye but extracting those into line features is proving difficult so i was wondering if anyone knows a good algorithm or method that can handle the large images that Landsat provides? I've tried Grass's r.thin but this doesn't seem to work.


texts: QGIS plugin : connect to polygonselect Maptool signals

texts:
Working on a new plugin I am trying to connect a custom function to the mouse right-clic signal while using the polygon selecting maptool.
I have activated the maptool through :
def run(self):
    self.iface.actionSelectPolygon().trigger()

And would now like to connect my selectSite function to a canvasReleaseEvent or similar signal. From an interesting tutorial here I adapted the code and tried the line below which failed :
def run(self):
    self.iface.actionSelectPolygon().trigger()
    result = QObject.connect(self.canvas.mapTool(), SIGNAL("canvasReleaseEvent()"), self.selectSite)

def selectSite(self):
    QMessageBox.information( self.iface.mainWindow(),"Info", "in selectSite function" )

Is there another way around by connecting to a more generic Qt mouse click signal ?
Any ideas or hints on how I can achieve my goal ?


texts: How to connect particular sets of points with (poly)lines?

texts:
I have a timestamped point-table representing GPS-tracked places of different cars. To distinguish the cars there is a car-id for every point. So, a simplified version is looking like that:

Now, I want to write a query for connecting the places of each car chronologically with lines so that the track of each car can be displayed. The table is sorted chronologically, so it can be worked with the id's instead of the timestamps.

How can I solve it?


texts: CartoDB basemaps using CartoDB.js with createLayer

texts:
I've created a map on the CartoDB site and now wish to add some interactivity using CartoDB.js.

I'd like the interactive map to have the same basemap as on the CartoDB site, Nokia Hybrid Day.

The createVis method automatically loads the CartoDB basemap. But I need to use the createLayer method, which does not.

Could you please tell me how I can load the basemap I'm using on the CartoDB site for use with the createLayer method?


texts:
I am using the QGIS (3.0.3) Advanced Digitizing Panel. I would like to control the spaping characteristics. The QGIS documentation states that by selecting the "gear" drop down icon...

  You can also control the snapping to features. The options are:

  • Do not snap to vertices or segments
  • Snap according to project configuration
  • Snap to all layers

However, when I select the "gear" the drop down menu does not contain these options. How do I control the snapping features in the Advanced Digitizing Panel in QGIS 3.0.3?


texts: Earth engine: function "is not defined in this scope"

texts:
I am trying to use in Google Earth Engine the function area() or centroid() on a simple polygon. According to the documentation, there seems to exist a method area() for objects ee.Geometry.Polygon. However, I get an error message with this simple code:

var polygon = ee.Geometry.Polygon([
  [[-5, 40], [65, 40], [65, 60], [-5, 60], [-5, 60]]
]);

var a = area(polygon);
var aa = polygon.area();

Gives message:

  "area" is not defined in this scope.

What am I doing wrong?