texts: Transform Google Earth Engine script to Python with Landsat 8 temporal data

texts:
I am doing a study for my work and I need to adapt my Google Earth Engine script in Python, but I am encountering some problems.
Here is a Google Earth Engine script:
 //Importing image and geometry:
var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
//    geometry = /* color: #d63000 */ee.Geometry.Polygon(
//        [[[-80.92489760146748, 25.433457120928352],
//          [-80.64474623427998, 25.488013471687964],
//          [-80.57882826552998, 25.710940372707608],
//          [-81.02377455459248, 25.770317250349557],
//          [-80.95236342177998, 25.552457242621447]]]);

//Filtering date, polygon, and cloudiness
    var image = l8.filterDate ('2010-09-01', '2021-12-31')
                .filterBounds (geometry)
                .filterMetadata ('CLOUD_COVER', 'less_than', 1);

    //NDVI calculation:
    var ndvi_func = function (i) {
      var ndvi = i.normalizedDifference (['B5', 'B4']).rename ('NDVI')
      return i.addBands(ndvi);
    }

    var image_ndvi = image.map(ndvi_func);

    //Calculating year wise NDVI
    var year = ee.List.sequence(2010,2021);
    var year_func = function(y){
      var range = ee.Filter.calendarRange (y, y, 'year');
      return image_ndvi.select('NDVI').filter(range).mean().set ('Year', y)
    };
    var yearwise_ndvi = ee.ImageCollection(year.map(year_func));
    print (yearwise_ndvi);
    Map.addLayer (yearwise_ndvi)

//Creating time-series chart:
var chart = ui.Chart.image.series ({
  imageCollection: image_ndvi.select('NDVI'),
  region: geometry,
  reducer: ee.Reducer.mean(),
  scale: 30})

print(chart);

Figure:

The script above presents an NDVI time series for a region and I need to do it in Python. Here is a script that shows the errors in Python:
l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
geometry = ee.Geometry.Polygon([[[-80.92489760146748, 25.433457120928352],
                                 [-80.64474623427998, 25.488013471687964],
                                 [-80.57882826552998, 25.710940372707608],
                                 [-81.02377455459248, 25.770317250349557],
                                 [-80.95236342177998, 25.552457242621447]]]);

#Filtering date, polygon, and cloudiness
image = l8.filterDate ('2010-09-01', '2021-12-31')
image  = image.filterBounds (geometry)
image  = image.filterMetadata ('CLOUD_COVER', 'less_than', 1);

#NDVI calculation:
def ndvi_func(i):
    ndvi = i.normalizedDifference (['B5', 'B4']).rename ('NDVI')
    return i.addBands(ndvi)

image_ndvi = ndvi_func(image.map)

#Calculating year wise NDVI
year = ee.List.sequence(2010,2021);
def year_func(y):
    range = ee.Filter.calendarRange (y, y, 'year');
    return image_ndvi.select('NDVI').filter(range).mean().set ('Year', y)

yearwise_ndvi = ee.ImageCollection(year.map(year_func));
    print (yearwise_ndvi);
    Map.addLayer (yearwise_ndvi)

#Creating time-series chart:
chart = ui.Chart.image.series ({imageCollection: image_ndvi.select('NDVI'),
                                region: geometry,
                                reducer: ee.Reducer.mean(),
                                scale: 30})

print(chart)

The problem is that the functions do not load the data.
How can I load NDVI data by date into a DataFrame?
I would also like to have daily NDVI instead of monthly / yearly, but I'm not getting it on Google Earth Engine.


texts: Save/Push drawn features in postgreSQL(PostGIS) using OpenLayers and GeoServer

texts:
I am working on a web-based GIS application using OpenLayers, Geoserver and PostgreSQL (PostGIS). Till now I am able to draw a polygon on the front end but drawn shape does not store in Database PostgreSQL (PostGIS). As I refresh the webpage, the drawn feature disappear because the ajax query is not properly working.  Need expert's input to store the data in DB PostgreSQL (PostGIS). Below is the code:
var formatWFS = new ol.format.WFS();
var formatGML = new ol.format.GML({
    featureNS: 'http://localhost:8080/geoserver/DBdata/wfs',
    featureType: 'wfs_geom',
    srsName: 'EPSG:3857
});

var s = new XMLSerializer();
// read data from the server
var sourceWFS = new ol.source.Vector({
    loader: function (extent) {
        $.ajax('http://localhost:8080/geoserver/DBdata/wfs', {
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typename: 'DBdata:district',
                srsname: 'EPSG:3857',
                bbox: extent.join(',') + ',EPSG:3857'
            }
        }).done(function (response) {
            sourceWFS.addFeatures(formatWFS.readFeatures(response));
        });
    },
    strategy: ol.loadingstrategy.bbox,
    projection: 'EPSG:3857'
});

var layerWFS = new ol.layer.Vector({
    source: sourceWFS
});

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
          source: new ol.source.OSM()
         }),
        layerWFS
    ],
    view: new ol.View({
        center: ol.proj.fromLonLat([71,30]),
        zoom: 4
    })
});

var interaction = new ol.interaction.Draw({
    type: 'Polygon',
    source: layerWFS.getSource()
});

map.addInteraction(interaction);
// post query
interaction.on('drawend', function (e) {
    $.ajax('http://localhost:8080/geoserver/DBdata/wfs', {
        type: 'POST',
        dataType: 'xml',
        contentType: 'text/xml',
        data: s.serializeToString(formatWFS.writeTransaction([e.feature], null, null, formatGML))
    }).done();
});

As of my understanding, post query should be executed and insert the drawn feature in the database but it is not happening.


texts: Arcgis javscript calculating lines in a place

texts: I have two mapservice. Fist service includes parcels and second includes water pipe lines. I want to calculate total pipelines length in a parcel. How can I calculate this via arcgis javascript API. For example 2201 number parcel includes 284 meters.


texts: Delimited text layer not showing up on map

texts:
I am using QGIS 3.14.15 on MacOS. I have added a delimited text layer (csv file converted from xlsx) as a new shapefile. It shows up on the "Layers" section at the bottom left but not on the map. If I click "Zoom to layer," it shows up but everything else disappears. When I undo it, the layer disappears again. Made sure coordinates were represented on X and Y axis. Any ideas?
adding delimited text layer
[


texts:
I am trying to reproject some climate data that I downloaded from http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/dcpInterface.html#Welcomein

The BCSD-CMIP5 data is found under the tab "Projections: Complete Archives".

The files are in netCDF format, but am running into some issues.  The original raster data (let's call it rast1)  when imported has the extent

xmin 235.25
xmax 293
ymin 25.125
ymax 52.875

and a crs of +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

The other objects (both sp and raster) I am working with contain negative values for xmin and xmax

xmin        : -122.8021
xmax        : -75.38249
ymin        : 31.78881
ymax        : 47.15952

When I try to reproject rast1 using the other objects (rast2), the extent does not change and even though the reprojected rast1 is in the same crs, the extents don't line up.

rast_new = projectRaster(rast1, rast2)

So I had to do the following:

rast_new = projectExtent(rast1,rast2)
rast_new = projectRaster(rast1,rast_new)
values(rast_new) = values(rast1) ### This is necessary because all the values were NA

However, even though rast_new is in the same projection as rast2

+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83
+units=m +no_defs +ellps=GRS80+towgs84=0,0,0

the plotting isn't the same. See image

Just wondering what I can do to correct this? Needless to say, I am relatively new to understanding projections in R ( and projections in general!)

I have exported the two rasters and viewed in ArcGIS, to see if it was a plot() issue in R, but the images are the same as above, yet the projections are the same (according to ArcGIS 10.2). Even when I reproject the "bad" raster in ArcGIS, the projection does not appear to be in Albers Equal Area. When I import directly from the netCDF file into ArcGIS, the projection looks correct.


texts: Load postgis database with PK from two fields in QGIS

texts:
I have a PostGIS database having a Primary Key from two fields id_cat and id_cdad. When I load this database in QGIS in my computer under Linux there are no problem. I unmark "Select ID" and I can mark more than one field in the combo box of Primary Key because there are a checkbox that let me do that. (QGIS 2.12 under Linux Mint 17.0)

Now I want to install my database in a computer under window, but there aren't the checkboxes close to the name of fields, so I can't to mark more that one field for define the primary key.
(QGIS 2.10 under windows 7)

Any suggestion?


texts: Make array_append in trigger with PostGIS intersection

texts:
I've 2 layers : linestring with column id_line and polygone, drawn by user.
I trying to get id_line on linestring who intersects polygon.
I've got it in qgis with :

aggregate(
layer:= 'aep_cana',
aggregate:='concatenate',
expression:=id_cana,
concatenator:='; ',
filter:=intersects($geometry, geometry(@parent)))

in default value, It's not really an array, more like a list but it's ok for me.

So I try to make the same thing when I create a polygon, automatically, in my database with a trigger (because in QGIS, users can do anything they want).

CREATE OR REPLACE
FUNCTION asst.transf_id_line() RETURNS TRIGGER LANGUAGE plpgsql AS $function$
BEGIN
    UPDATE asst.test_emprise SET
    list = array_append(list, (SELECT c.id_line FROM line."id_line" AS c WHERE st_within(NEW.geom, c.geom))) ;
END;
$function$
;

Doesn't work, I write this code just for idea.
But I try some simple array_append in trigger, but it makes loop...

So I find WITH RECURSIVE like in
Postgis: aggregate function for intersection
but I can't apply it to me.
It this a good way or not? How can I use array_append in a trigger?

After 2 hours, I can make

WITH RECURSIVE idsel AS (
SELECT
    te.id AS idi,
    ARRAY_AGG( c.id_line::TEXT) AS aa
FROM
    asst.test_emprise AS te
CROSS JOIN asst."line" AS c
WHERE
    ST_INTERSECTS(c.the_geom,te.geom)
GROUP BY
    te.id )
--SELECT * FROM idsel;
UPDATE asst.test_emprise AS op SET list = idsel.aa FROM idsel WHERE idsel.idi = op.id;

But I can't transform in trigger.


texts: Field calculator script to work out first 5 dominant area fields?

texts:
I have 5 fields of different landcover types with areas of each and found below:

> PRC_Water| PRC_Bare| PRC_Urban| PRC_LVeg| PRC_LVeg|
>     5        20         40        20         15

I have 5 more fields of dominance [Dom1, Dom2, Dom3, Dom4, Dom5] which I would like to populate with the field_name of the landcover type. I would like to know if it is possible to write a script/conditional statement in the field calculator to return the field name of the '#th' largest area, and if there is a way to return both names if the areas are equal??

I have limited knowledge of python, but am puzzled how I can rank these areas by returning a field name?


texts: Creating a CoverageView with GeoServer API

texts:
Using GeoServer Rest API, I have uploaded a 4 band Geotiff image in a coverage store. But I can't seem to find the right request to create a CoverageView from this geotiff, choosing which band this CoverageView will have, although I can do it from GeoServer interface. Could anyone point me in the right direction ?
Like in this post, I created a coverage store and uploaded my geotiff.
To create the coverage view, I send this query, using httr package in R :
r = POST(
     url = "http://localhost:8080/geoserver/rest/workspaces/my_workspace/coveragestores/my_coverage_store/coverages",
     config = authenticate(id, password),
     add_headers(.headers = c("Content-Type" = "text/xml")),
     body = body)

and the xml body :
<coverage>
    <name>my_coverage_store</name>
    <nativeName>my_coverage_store</nativeName>
    <metadata>
        <entry key=\"COVERAGE_VIEW\">
        </entry>
    </metadata>
</coverage>

And with content(r) I get this error :
"java.lang.String cannot be cast to org.geoserver.catalog.CoverageView"

Alternatively, another answer for my problem would be to define a sld or css style where I can choose which band to display and specify that no alpha band should be used. So far I wasn't able to do that since in the openlayers preview of my layer using JPEG is ok, but with PNG the 4th band is set as an alpha band.
Thanks !


texts: How to Get a Vector Layer to add to a QGIS Project

texts:
I used QGIS a little bit a few years ago and used some layers for my city that I had. I have a simple project that I would like to use it for and am having a little problem. I would like to simply get an orthophoto of a parcel that we recently purchased. I found this thread.

It is practically what I want to do.  I just want the data for another township in the same state.

I can set up the coordinate system and know the WMS layer I would like to use for the orthos, but, the instructions say to add a vector layer of the specific area I'm interested in.  The example gives Acme Township as an example.  The instructions don't tell me where I can find the layer so I can add it.  In my case I am looking for Mecosta Township, MI.

Can anyone tell me where I might be able to find some sort of shapefile for Mecosta Township. If the layer includes parcels that would be ideal.


texts: QGIS (2.6.1) Processing Toolbox InitGui error after changing /Users folder on Mac OSX

texts:
Using: Mac OSX 10.9.4
       QGIS 2.6.1

I had previously installed QGIS and had everything running fine. My home directory was /Users/rowanbado, but for another application I had to add an intermediate directory such that my home directory is now stored in /Users/bado/rowan. Now when I open QGIS I get the following error:

Couldn't load plugin processing due an error when calling its initGui() method

Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 217, in startPlugin
plugins[packageName].initGui()
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/ProcessingPlugin.py", line 54, in initGui
Processing.initialize()
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/core/Processing.py", line 145, in initialize
Processing.modeler.initializeSettings()
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/modeler/ModelerAlgorithmProvider.py", line 55, in initializeSettings
ModelerUtils.modelsFolder()))
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/modeler/ModelerUtils.py", line 46, in modelsFolder
mkdir(folder)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 127, in mkdir
mkdir(head)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 127, in mkdir
mkdir(head)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 127, in mkdir
mkdir(head)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 129, in mkdir
    os.mkdir(newdir)
OSError: [Errno 13] Permission denied: '/Users/rowanbado'

Which is particularly confusing because that directory no longer exists. I have tried the following workarounds:

-Uninstall and reinstall QGIS (as well as all frameworks, cairo,gdal,etc, and delete ~/.qgis2 folder). Didn't do anything.

-Comment out the line in 'ModelerUtils' (Line 46) which appears to be the source of error. I no longer receive an error message from 'ModelerUtils.py' (so I assume the module is running fine) but instead receive the same message as above except with SystemUtils.py, which I assume is the next module to be initialized. Clearly the issue is with the system.py module.

-Create a symbolic link in /Users/bado/rowan to /Users/rowanbado. Now everything works fine and I have no problems, but this is sort of a 'dumb' work around in that it treats the symptom and no the problem.

It seems that the Processing toolbox has stored a line with my old home directory path in some internal file, and now fails upon initialize and therefore cannot change this line as I assume it should. When I start up the (now working with a symbolic link) toolbox, then remove the symbolic link, I am right back to the same error.

Any developers out there have some idea what's going on and how I can address this? Seems to be a trivial issue to someone who understands the system very well.


texts: Trying to fill polygon fills the entire map

texts:
I have been given some shp files to be displayed using Mapserver. Until now I was able to find the correct projections, display the map, and show labels. However, I have problems with polygon colors. When I set COLOR to any value, the entire map is filled with that color.

Is it a problem with the shp file, or what am I doing wrong?

For example:

LAYER
    NAME mylayer
    TYPE POLYGON
    DATA "mylayer.shp"
    STATUS DEFAULT
    CLASS
      NAME mylayer

      STYLE
        COLOR 255 0 0
        OUTLINECOLOR 0 0 0
      END
    END
  END

This will fill the entire map with red color, polygons have black outline.


texts: Calculating line length within mulitple overlapping polygons using ArcGIS Desktop?

texts:
My question is similar to Calculating edge length of polygon within overlapping polygon using ArcGIS Desktop?

However, my issue differs in that I have multiple overlapping polygons (upwards of 7) that share the same line.

Here's an example:

Polygons represent 50 meter buffers around a point and the line represents the edge of another polygon that I've clipped.

I'm trying to calculate the length of the blue line contained within multiple (yellow) polygons. The issue I'm running into is that I'm not able to calculate the length of the line within each individual polygon. I've tried Spatial Join and Intersect which haven't worked. Spatial Join results in the total length of the line regardless of polygon. Intersect segments the line by polygon border, therefore cutting off portions that should be included within a given polygon, i.e. length isn't calculated for a polygon that is overlapped.

I am using ArcGIS Desktop 10.6.


texts:
Is it possible to define point grid spacing by attribute in QGIS?

EDIT: I really dunno is this even possible. Added a picture where I try to demonstrate the situation.

So the problem is that I want to generate a grid where I can create a grid that takes the spacing value from an attribute table like the one below. And also that the grid is only placed inside the polygons. Normally I do this so that i create a grid and the I select the points that are located inside a specific polygon and delete all the points that are not inside the wanted polygons.

(I'm not very familiar to using the python console,  but every answer is welcome.)

Example:

(Simple question, but I didn't find any answers to it.)


texts: Add properties name label or id to a feature

texts:
I'm new to OpenLayers 3 and I'm trying to add a property (id, name or label) to a feature to get it back on click later, to identify that particular feature.

Is it possible?

I haven't found any answer doing this yet. Something like:

var point_feature1 = new ol.Feature({ });
point_feature1.setProperties('description', "description 1");
point_feature1.setId(666);
...

And then onclick:

map.on('singleclick', function(evt) {
  var feature = map.forEachFeatureAtPixel(evt.pixel,
    function(feature, layer) {
      console.log("feature",feature.getId());
      return [feature, layer];
    });
});


texts:
How can I find the closest feature(coordinates) in a json file from a set lat lon using openlayers.
Currently, I get the centroid point(lat,lon) from the selected polygon feature ok, but I  Would also like to return the closest feature "sid" value from the json file below.

Here is the contents of the json file:

 { "type": "FeatureCollection",
    "features": [
      {
  "type":"Feature",
    "geometry":{"type":"Point","coordinates":[146.9509,-36.069]},
      "properties":{"name":"ALBURY AIRPORT AWS","sid":72160,"rain":[0,0,0,0,0,0,0,1.2,0,0]}
},
{
  "type":"Feature",
    "geometry":{"type":"Point","coordinates":[142.0867,-34.2358]},
      "properties":{"name":"MILDURA AIRPORT","sid":76031,"rain":[0,0,0,0,0,0,0,0,0,0]}
},
{
  "type":"Feature",
    "geometry":{"type":"Point","coordinates":[142.3158,-35.0694]},
      "properties":{"name":"OUYEN (POST OFFICE)","sid":76047,"rain":[0,0,0,0,0,0,0,0,0,0]}
},
{
    "type":"Feature",
      "geometry":{"type":"Point","coordinates":[142.004,-35.1201]},
          "properties":{"name":"WALPEUP RESEARCH","sid":76064,"rain":[0,0,0,0,0,0,0,0,0,0]}
},
{
  "type":"Feature",
    "geometry":{"type":"Point","coordinates":[143.5416,-35.3766]},
      "properties":{"name":"SWAN HILL AERODROME","sid":77094,"rain":[0,0,0,0,0,0,0,0,0,0]}
}]}


texts: Why would running LocateFeautresAlongRoutes GP tool return null?

texts:
I'm currently using the LocateFeaturesAlongRoutes tool wrapper for C#/.NET for a custom tool I'm developing for ArcMap 9.3.1.  I've been able to run this tool in ArcMap outside of code in a tool with the same files as parameters.

Unfortunately when I put it in code I keep getting null for the IGeoProcessorResult.  Is there any reason why I would get a null object for the result instead of just saying it failed? Also, is there anything wrong with the below code?

The relevant block of code is below.

LocateFeaturesAlongRoutes lfar = new LocateFeaturesAlongRoutes();
lfar.in_features = @"C:\GIS_Temp\CONGPROP.mdb\CONGPROP2";
lfar.in_routes = @"C:\GIS_Temp\Routes.mdb\routes";
lfar.route_id_field = "ROUTE_ID";
lfar.out_table = "eventTable.dbf";
lfar.out_event_properties = "ROUTE_ID LINE FMEAS TMEAS";
lfar.in_fields = "NO_FIELDS";

Geoprocessor geoprocessor = new Geoprocessor();
geoprocessor.TemporaryMapLayers = false;
geoprocessor.OverwriteOutput = true;
geoprocessor.AddOutputsToMap = true;
IGeoProcessorResult result = (IGeoProcessorResult) geoprocessor.Execute(lfar, null);

Any suggestions?


texts: ogr2ogr import fails with "ERROR 1: GeoJSON object too complex"

texts:
I'm trying to import a large (8.8GB) GeoJSON file to a Postgres 9.6 database with the following command:

ogr2ogr -f "PostgreSQL" -nln final_diffs PG:"host=localhost dbname=inspire user=me" diffs.geojson

But I get lots of warnings like this:

ERROR 1: GeoJSON object too complex
ERROR 1: GeoJSON object too complex
More than 1000 errors or warnings have been reported. No more will be reported from now.

And the process runs without errors, but at the end, nothing has been imported into the Postgres database.

Any ideas? I've not seen this warning before, and the only online reference I can find is in the GDAL source code.

All I get from ogrinfo is the above warnings, then:

Layer name: diffs
Geometry: Unknown (any)
Feature Count: 0
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]


texts: Amending existing plugins to match Qgis 2.16+ fTools removal

texts:
Part-time development started about a year ago to create an add-on for Qgis.  Many of this plugins capabilities are written with fTools commands.

I am wondering if there is a recommended document to help developers migrate syntax to the new processing functions included with 2.16+.

Although fTools is still available through the Plugin Manager, rewriting the plugin is best for legacy with the Qgis programs future.


texts: Can I import symbology from an .axl file into ArcMap?

texts:
I'm setting up new Arcserver services, and some users would like the symbology in the new services to match the old IMS services.  With that in mind, is it possible to either:

a) import the symbology from an axl file into ArcMap
b) convert the axl file to an avl file so the symbology can be imported

Thanks in advance for any help, any ideas appreciated.


texts: Training courses in geostatistics and spatial modelling

texts:
Possible Duplicate:
  What are sources for current geo-statistical analysis classes/events: open or proprietary

Does anyone know of or have been on any training courses that specifically train in geostatistics or spatial modelling techniques with reference to the mapping (ArcGIS) world?  E.g. Kriging, interpolation etc..


texts: How to make a polyline with two-sided arrows show two colors

texts:
I am very new to ArcGis development and so bear with me if I asked a silly question. I am using ArcGis 10.2 on C#.
I have two features (A and B) on the map. I need to place two arrowed lines - one from A to B (pointing B), and one from B to A (pointing A). Both sets are expected to show different colors, which will convey different meanings to the users.
Currently, I have a PolylineClass for the line and ArrowMarkerSymbolClass for the arrow, which are eventually packed into LineElementClass before adding to one of the map layers. This is done twice - once from A to B and another time for B to A. The problem is that the two sets of arrowed lines are overlapping. Because the ability to see the colors of the two sets of arrowed lines is the most important function here, I need to change how it is implemented.
I have two ways in mind, but I have no idea if it is achievable.

Draw two separate arrowed lines side-by-side. I totally have no idea how this is possible because the two lines share the same two coordinate points on the map.
Draw half of each arrowed line to a color.

Current Implementation:

Method 1:

Method 2:

So, is there any way I can achieve either #1 or #2?
Update
Posting some codes so that it's easier to understand what I'm doing.
    private void DrawArrowedLine(Place from, Place to, RgbColor color)
    {
        ESRI.ArcGIS.Geometry.PolylineClass polyLine = new ESRI.ArcGIS.Geometry.PolylineClass();
        polyLine.SpatialReference = MapHelper.GetSpatialReference();

        ESRI.ArcGIS.Geometry.Point fromPoint = new ESRI.ArcGIS.Geometry.Point();
        fromPoint.SpatialReference = MapHelper.GetSpatialReference();
        fromPoint.X = from.Location.Longitude;
        fromPoint.Y = from.Location.Latitude;

        ESRI.ArcGIS.Geometry.Point toPoint = new ESRI.ArcGIS.Geometry.Point();
        toPoint.SpatialReference = MapHelper.GetSpatialReference();
        toPoint.X = to.Location.Longitude;
        toPoint.Y = to.Location.Latitude;

        polyLine.AddPoint(fromPoint);
        polyLine.AddPoint(toPoint);

        ESRI.ArcGIS.Display.ArrowMarkerSymbolClass arrow = new ESRI.ArcGIS.Display.ArrowMarkerSymbolClass();
        arrow.Color = color;
        arrow.Size = 12;
        arrow.Length = 12;
        arrow.Width = 8;

        ESRI.ArcGIS.Display.CartographicLineSymbolClass cartoLineSym = new ESRI.ArcGIS.Display.CartographicLineSymbolClass();
        cartoLineSym.Color = color;
        cartoLineSym.Width = 3;

        ESRI.ArcGIS.Display.SimpleLineDecorationElementClass lineDecElem = new ESRI.ArcGIS.Display.SimpleLineDecorationElementClass();
        lineDecElem.AddPosition(1);
        lineDecElem.MarkerSymbol = arrow;

        ESRI.ArcGIS.Display.LineDecorationClass lineDec = new ESRI.ArcGIS.Display.LineDecorationClass();
        lineDec.AddElement(lineDecElem);

        cartoLineSym.LineDecoration = lineDec;

        ESRI.ArcGIS.Carto.LineElementClass lineElem = new ESRI.ArcGIS.Carto.LineElementClass();
        lineElem.Symbol = cartoLineSym;
        lineElem.Geometry = polyLine;
        lineElem.Name = CommsLinkName;
        lineElem.SpatialReference = MapHelper.GetSpatialReference();

        // ArrowLayer is ESRI.ArcGIS.Carto.CompositeGraphicsLayerClass
        MapHelper.GetArrowLayer().AddElement(lineElem, 0);
    }


texts: QGIS Web Client 2 & QGIS Server features break after adding new project source

texts:
QUESTION:

Ubuntu 18

I have problems with QGIS web client 2 features after install & configure w/ QGIS server.

I must be doing something incorrect with the project html link if I am seeing these issues i detail below, because in the demo all works correctly and the html links do not include wms request.

I see the html links in the demo themesConfig.json are like this:

http://qwc2.sourcepole.ch/wms/themes/test

So no capabilities request or WMS?

When I leave out the wms request from the html yarn does not want to reload the theme even after symbolic links for project source, wms xml, & mapserv.fcgi

I can only get the map to appear when I do the html link like this:

http://qgis.demo/cgi-bin/projects/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

I assume this is the core of the issues I am seeing. What am I missing?

ISSUES:

1: When I use print function, I get 500 internal server error after pressing "print" button.

test project has a print composer profile setup.
something in themesConfig.json or Config.json missing or incorrect possibly

2: When I click on the map the pop up window that would show map feature details in the demo quickly disappears and I am sent to a blank white screen and I must refresh browser to get the map back.

this feature was working with demo project sources & broke when I added my source

3: When I export RASTER & DXF I get "unable to connect localhost:5000"

same in demo & after I added my project source.
I see in config.json i need to configure proxy, but details are limited.

4: When I access QWC2 from another system on the network, I am unable to see the test project. The thumbnail and project layers are not showing up.

DOCS USED:

Followed steps here for QGIS server install, sections 12.1 & 12.2:
https://docs.qgis.org/2.18/en/docs/training_manual/qgis_server/index.html
Followed steps here for QGIS web client 2 after 12.2.6 QGIS-server instructions said to use the older QWC version:
https://github.com/qgis/qwc2-demo-app/blob/master/doc/QWC2_Documentation.md

RESULTS:

QGIS install was successful
can send WMS get capabilities and map requests.

QWC2 demo install was successful
all features working with demo project sources (expect export RASTER & DXF, I get "unable to connect localhost:5000")

I modified the themesConfig.json & the test project (world.qgs) displays in QWC2 locally on server:

http://qgis.demo/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&map=/home/qgis/projects/world.qgs

or

This (w/ symbolic links):

http://qgis.demo/cgi-bin/projects/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities


texts: Specific projection won't work, while others do

texts:
When I do this in JavaScript, using OpenLayers 3, it works:

    ol.proj.transform([153.47, -28.11], 'EPSG:4326', 'EPSG:3857');

I realize that this is because OpenLayers already knows both projections and is able to transform between them.

I want to be able to project from 'EPSG:28356' (GDA94 zone 56) to 'EPSG:3857'. So, I've installed and imported proj4js, got the appropriate string from http://spatialreference.org/ref/epsg/gda94-mga-zone-56/proj4js/ and it gives me:

"+proj=utm +zone=56 +south +ellps=GRS80 +towgs84=200,0,0,0,0,0,0 +units=m +no_defs"

So, in my code:

    proj4.defs("EPSG:28356", "+proj=utm +zone=56 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");

    let projection = new ol.proj.Projection({
      code: this.props.item.data.code,
      units: this.props.item.data.units,
      extent: this.props.item.data.extent
    });

    ol.proj.addProjection(projection);

    let c1 = [543100.508297332, 6880073.9025688];
    let c2 = ol.proj.transform(c1, 'EPSG:28356', 'EPSG:3857');

But this doesn't work. It didn't surprise me to get the exact same coordinates (somewhat rounded) back in c2, since the string states '+towgs84=0,0,0,0,0,0,0'.

I can also get the correct result from https://epsg.io/transform#s_srs=28356&t_srs=3857&x=543100.5082973&y=6880073.9025688 which would be [17080772.24, -3274771.77].

Does spatialreference.org give me the wrong proj4 string for this projection and if so, where to get a correct one? Or am I making some other mistake in trying to get this projection to work?


texts: Update programmatically BBOX of WMS Geoserver layer

texts:
I create programmatically a layer in PostGIS and Geoserver. I give the option to the user to choose the geometry and the attributes and then I create the table in PostGIS and add it in an existing store in Geoserver.

The layer, as initialized, it gets a global Bounding Box.
Then the user is able to add specific data. What I want is to be able to recalculate the Bounding Box depending on the data, so when I show the map on Leaflet to show only the selected area.
I need to do that programmatically using the Rest API.

I haven't found anything like this for single layers (not layergroups) and also for WMS layers.

Does anyone has any experience on it?


texts: Assigning CRS using las2las from libLAS?

texts:
My las file had a spatial reference of "None." I tried to assign the CRS (in which data were collected) with the following command line:

las2las flightline_002.las --a_srs EPSG:32632
lasinfo --no-check flightline_002.las

Lasinfo still reports it as "None." See screenshot below, line 'Spatial Reference':

Those extent coordinates are UTM meters.

I am working through the examples on http://www.liblas.org/start.html#osx, with libLAS 1.7.0, GeoTIFF 1.4.0 GDAL 1.11.1, Mac OS 10.9.4. Any thoughts?


texts: How can I figure out whether a sphere is intersecting the Earth between a given lon/lat range?

texts:
Given a lon/lat range that represents a tile on the Earth's surface, I need a computationally efficient way to figure out if a spherical volume defined by a point (assume ECEF coordinates) and a radius intersect with the tile.

It can be assumed that the maximum tile side length won't exceed 22.5 degrees  and that all lon/lat area coordinates are aligned along traditional map tiling schemes (tiles won't cross the antemeridian, etc).

In the following picture, the blue ellipsoid represents Earth and the green sphere represents a volume of interest.

Given a lon/lat area, you can generate ECEF coordinates for the 4 extreme points and check whether any of these points lie within the sphere through a simple distance calculation. However, this method doesn't always work since you're ignoring curvature. The following picture shows this, with the tile shown in black and a possible sphere of interest in green.

You can add more points to test with the lon/lat tile (ie. check the center of the tile along with the 4 corners) but you're just trying to improve an approximation, and this will fail based on the position and size of the sphere, so I want a more analytical solution. I'm looking for a solution I can implement myself in code and not for an application or library.

I'd like to use an ellipsoid representation of Earth if possible, but if that ends up being too complex, representing the Earth as a sphere is okay too.


texts: How do I merge a clipped raster on top of another raster is QGIS?

texts:
I am trying to merge two rasters (geotiffs) in QGIS.

One raster is the original & the other is clipped to a specific region of the original raster and 0 applied to null area to make transparent. See pictures below.

Original:

clipped:

Both together, clipped on top:

I want to merge the clipped raster on top of the original raster to make one file where the original raster is visible in the null (transparent) area of the clipped raster.

Depending on the order I list the files in the gdal command generated by merge function in QGIS I get a merge result where you cannot see the clipped raster or I get a merge result where the original raster is not visible in the null area of the clipped raster.

What do I need to do to have the clipped raster merge on top of the original raster and be able to see the original raster in the null area of the clipped raster?

Some of the questions I have looked through that seemed related:

Which value does raster merge chose in overlapping area?

How to merge two or more rasters which overlap?

Merging raster in QGIS?

Per request of 1st responder here is the gdal merge command generated in qgis:

gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /home/qgis/Desktop/merge1.tif /home/qgis/Desktop/41_36_1.tif /home/qgis/Desktop/41_36_2.tif

I have tried with the files in both orders.

Here's a link to download the files I am working on.

https://drive.google.com/open?id=17wyRKbkDE1twWBSMDsY2upY7cuoHmOq1


texts: Selecting Points by their Spatial Location to an Attribute of a Layer

texts:
I have a large point layer of addresses within a district. Each address point has a field called "Street External Key [SEK]" which is tied (meaning identical to) the corresponding "Street External Key [SEK]" field that I have in another layer, a street layer for the district.

The idea is to be able to associate each address to the street off of which it is accessed (ex. 111 Main Street will have a SEK value of 22222, which is the SEK value of the street that it is right next to).

However, some errors have popped up where an address was tied to a SEK of a street miles away.

How would I be able to perform a search to get all addresses where the SEK is from a street that is, say more than 250 ft away, which means its probably an error?

I essentially want to select all points that are 250 ft or greater from a street (which to my understanding is select by location), where the street is defined by an attribute [SEK] that both features have values for (hence the attribute qualifier i need in addition to the select by location).

I looked at this Q&A, but have been unable to make the connection to my situation:

Using Select By Attribute and Select By Location simultaneously?


texts: How do you merge multiple polyline features into one contiguous feature in QGIS?

texts:
I created a shapefile from public GIS data, that I wanted to display in Google Maps for presentation purposes. The original .shp used small line segments for roads. Maybe a picture will make sense:

As you can see, the feature class is made up on 164 little line segments. I want it to be one contiguous feature, so that when I import into Google Maps, it doesn't list every separate line segment. How can I do that?

For the record, I'm using QGIS 3.6.1 with GRASS 7.6.1.


texts: access attributes of one class for use in another class in python add-in

texts:
I want to get combobox attribute 'dem' and use it for further processing in class ExtentData. But I am not sure how to get attributes from one class to another for processing in pythonaddin. In following code, I want attribute 'dem' from class ComboBoxFrame right after if statement in class ExtentData. Any suggestion would be appreciative.

import arcpy
import pythonaddins
import wx

class ExtentData(object):
"""Implementation for ExtentData_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
        self.cursor = 5
        self.shape = 'Rectangle'

    def onRectangle(self, rectangle_geometry):
    """"""

    dlg = None
    extent = rectangle_geometry

    if extent != None:
        self.dlg = ComboBoxFrame()

class ComboBoxFrame(wx.Frame):
    def __init__(self):
        wx.Frame.__init__(self, None, -1, 'Input Data Selection', size=(350, 300))
        self.Bind(wx.EVT_CLOSE, self.OnClose)
        panel = wx.Panel(self, -1)

        DemList = ['DEM1', 'DEM2', 'DEM3']
        self.Dem = wx.ComboBox(panel, -1, "", (15, 30), wx.DefaultSize,DemList, wx.CB_DROPDOWN)
        self.Bind(wx.EVT_COMBOBOX, self.OnSelect, self.Dem)

        self.btnApply = wx.Button(panel, label="Apply", pos=(65, 200))
        self.Bind(wx.EVT_BUTTON, self.OnClose, id=self.btnApply.GetId())
        self.Show(True)

    def OnClose(self, event):
        self.Show(False)

    def OnSelect(self, event):
        dem = str(event.GetValue())

app = wx.PySimpleApp()
app.MainLoop()


texts: How to project pixel coordinates in leaflet?

texts:
I am trying to create a map based on an image which is 16384x16384 px, but I would also like to add markers at specific locations using pixel coordinates of this given image.

I created a tile layer, a map element and set maximum boundaries, so that I can't scroll out of the image, using this code:

var map = L.map('map', {
    maxZoom: 6,
    minZoom: 2,
    crs: L.CRS.Simple
}).setView([0, 0], 2);
var southWest = map.unproject([0,16384], map.getMaxZoom());
var northEast = map.unproject([16384,0], map.getMaxZoom());
map.setMaxBounds(new L.LatLngBounds(southWest, northEast));

L.tileLayer('tiles/{z}/{x}/{y}.png', {
    tms: true
}).addTo(map);

The points i'd like to add to the map are stored in an external file as GeoJSON, and they look like this http://www.de-egge.de/maps/terranigma/terraPoints.js

I load them using this snippet of code:

var terraPoints = new L.geoJson(points, {
    onEachFeature: function (feature, layer) {
        layer.bindPopup(
            "<b>Ort: </b>" + feature.properties["points[, 1]"]
        );
    }
});

map.addLayer(terraPoints);

But of course, they don't show up, because the reference systems don't match. The points use pixel coordinates, but my map uses geographic coordinates.

Is there a way to "unproject" the points while loading them?


texts: memory issue when trying to buffer/union large dataset using postgis

texts:
I am trying to generated a dissolved buffer of a large point dataset (Ideally 29 million points - Address data in Great-Britain, but I receive the data by chunks of 1 million points, so these could be used instead).

After trying (and obviously failing!) different methods using ArcGIS, I came across this article:
http://blog.cleverelephant.ca/2009/01/must-faster-unions-in-postgis-14.html

Since I have been willing to try postgis for a while a gave this a go. Obviously being an absolute beginner I struggled a bit to get the whole thing up and running, but I finally managed to import data and buffer it. I am struggling however with the ST_Union call, because of memory issues (I am using postgis 2.0.1 running on pretty old 32bits windows pc with little scope to change that in the near future).

I tried changing the memory settings of postgres, as described in the question below:
Queries returning very big datasets in PostGIS

I believe I also created a spatial index properly on my geometry column. As far as I can tell, it didn't help.

I managed to run ST_Union on a 500 000 points subset, which is already pretty big. I think at that point the best for me might be to try using the Declare/Fetch method mentioned in the above question, however being a complete beginner I am having difficulties.

I believe it is in my interest to keep the slices as big as possible (because the ST_Union algorithm is fast and optimised). I think I would therefore just need to find a way to "slice" my ST_Union() queries to output one multipolygon for each 500000 rows. I would then be just one ST_Union on the resulting multipolygons away from the holy grail...

It is also possible that my query is sub-optimal. Or, there might by another great simple solution that I just haven't thought about or found online yet...

Here is a query I am currently using:

    drop table if exists ab_buffer_unioned;
    create table ab_buffer_unioned (
    gid serial primary key
    , the_geom geometry
    );

    insert into
        ab_buffer_unioned (the_geom)
    select
        st_Union(the_geom)
    from
        (select the_geom from ab_buffer limit 500000) as absubset;

Does anybody have any idea on how to fix my issues, or modify that query to add a loop or a cursor?

Many thanks

UPDATE

Many thanks to Paul for his suggestion below. The GeoHash trick worked great as can be seen below (see how the cluster size varies with the population density!)


texts:
I am interested in a DTM made available via a WMS service. I do know how to add the DTM as layer in ArcGIS (10.1 advanced lisence), but I do not know if it is possible to save the DTM (or part thereof) locally as a raster dataset.

Is there any procedure to get that accomplished?


texts:
We’re working on an evaluation of parcel conformity within zoning districts for a city and need to calculate frontage for all parcels. We’re able to do this for all parcels with the exception of corner parcels because of there being frontages for two streets. We’ve tried several approaches to connect the parcel address with the correct frontage (facing the address’ street) without result.

Does anyone have experience with addressing this issue?


texts: Import ESRI ASCII raster with python to GRASS

texts:
I had to get used to GRASS the last two days because I can't set up my python IDE to work with my osgeo4w python and I find it incredibly clunky. I have a simple ascii raster - I think its origin is ESRI ArcGIS given the following header structure:

ncols         389
nrows         265
xllcorner     665140
yllcorner     4823160
cellsize      20
NODATA_value  -9999

All I want to do is a simple import with the GRASS python editor to perform simple raster analyes. None of the 'Import raster data' menu entries worked out (probably because I can't flag the input type?), nor with the python editor. I tried to use r.in.gdal:

import grass.script as grass
grass.run_command("r.in.gdal",
                input="AAIGrid_raster.asc",
                output="GRASS_raster.asc",
                flags="f")

But I don't know how to specify the import format (most likely AAIGrid), as the description here doesn't help me sufficiently (btw. I think there are errors too on that page). In case I'd use -f AAIGrid in the shell, how does that look in pythonian?


texts: How to Set a feature layer Editable? (without Editor widget)

texts:
I'm new to arcgis JS api, now I want to update the whole records in a specific field of the feature layer, but I don't want to use the Editor Widget, is it possible? How to set the layer editable?
I did some research about using ApplyEdit() method, but still confuse about it.
Can anyone give me some hints or code examples? Thanks in advance!


texts: conditional st_join and st_intersection in R (sf)

texts:
I just started working with the R sf package and am trying to rewrite a simple procedure from PostGIS to extract shared borders between polygons. In PostGIS, I usually do a "conditional" self-join that filters the rows before doing the intersection (e.g. to avoid intersecting a lineshape with itself):

 WITH table1 AS
 (SELECT
    id,
    ST_Boundary(polygon) as geom
  FROM table0)
 SELECT
   t1.id,
   t2.id,
   ST_Intersection(t1.geom, t2.geom) as geom
   FROM table1 t1
     LEFT JOIN table1 t2
     ON ST_Intersects(t1.geom, t2.geom)
     AND t1.id<t2.id

In R, I can get the same result by doing the intersection first and subsetting the results after:

library(sf)
library(raster)

## Load example data (raster package)
fr0 <- getData('GADM', country="FRA", level=0)
fr1 <- getData('GADM', country="FRA", level=1)

## Get border intersections between units
fr1 <- st_boundary(st_as_sf(fr1))
fr1.is <- st_intersection(fr1, fr1)

## Remove self-intersections
fr1.is <- fr1.is[fr1.is$ID_1<fr1.is$ID_1.1,]

Result:

This works fine, but it would be much more efficient to do a conditional join first and then do the intersection, as is usually done in PostGIS. Unfortunately, my R skills are not that great yet. Any suggestions on how I can do this?


texts: How to select distinct values from a column in Mapbasic?

texts:
I have a table in a mapinfo .tab file. Is there any way to select distinct values from the column?

I tried :-

SELECT distinct (some_column) from table_name into temp_table

But this query would select all values. I wanted to select only distinct value from that column. Is there any way to achieve it?


texts:
My question is about availabe QGis Plugins that could help me with the following issue:

I have a PostGis database that has one shapefile as a main table and a lot of tables containing the shapefile's metadata linked to it. I want to query the main shapefile using the information in all the different tables. The result of the multiple-table-query will be a selection on the shapefile / a new shapefile that contains the selection. I can do this by constructing a complex SQL Query in DB Manager and loading the selection as a new layer.

BUT, the future users of my database can't write SQL queries, that's why I am searching for a graphical interface that enables them to construct multiple table queries.

I wonder if I missed any possibility of querying PostGis databases in QGis with help of a graphic query builder. I know that there are several Plugins as:

DB Manager, that is maybe the most complete but needs to write SQL queries
PostGis Query Editor, that has the same problem: users need to write SQL
PostGis Query Builder, that offers a graphic query builder and helps the future user of my database but it is limited to spatial queries
I have read that there was the "RT Sql Layer Plugin" but it was abandoned when DB Manager was developed, it was the kind of query builder interface that I am searching for but I can't find the code online any more


texts: QGIS 1.8 crashes when exporting from print composer to image on Windows 7

texts:
My question is similar to QGIS Crashes when writing to PDF.
Today, when I try to export to image, QGIS crashes. I've tried opening a new composer to no avail.

This is the problem signature:
  Problem Event Name:   APPCRASH
  Application Name: qgis.exe
  Application Version:  0.0.0.0
  Application Timestamp:    4fdce845
  Fault Module Name:    QtGui4.dll
  Fault Module Version: 4.7.1.0
  Fault Module Timestamp:   4d5308cc
  Exception Code:   c0000005
  Exception Offset: 00053de0
  OS Version:   6.1.7601.2.1.0.768.3
  Locale ID:    1033
  Additional Information 1: 0a9e
  Additional Information 2: 0a9e372d3b4ad19135b953a78882e789
  Additional Information 3: 0a9e
  Additional Information 4: 0a9e372d3b4ad19135b953a78882e789

Is there something wrong with my QtGui4.dll? Do I not have enough memory allocated?


texts: Error connecting to specific WFS layer with QGIS

texts:
I am trying to connect to the following WFS using QGIS:

  http://mapservices.prorail.nl/arcgis/rest/services/Basis_topografie_002/FeatureServer

However this gives me the following QGIS error:

  Capabilities document is not valid

  error occurred while parsing reference

I am using QGIS version 2.14

Is my URL wrong? I do get information if I open this in my browser.


texts: Performing identity operation in PostGIS?

texts:
I would like to perform spatial identity operation (ArcGIS example here) using PostGIS. I have two polygon tables, one of which is the input layer and the other table is the overlay layer. The output should be the original input layer but broken up to pieces where the overlay layer intersects it. This could get complex as multiple overlay polygons can intersect the input polygons. The output table need to have all the original fields of the input table plus a column from the overlay table. This extra column will have the data from the overlay table where they intersect and a fixed value like -1 or NA in non-overlapping areas.

How can this be done in PostGIS?

I am looking for an efficient solution as my input table has about 400k and my overlay table has about 20 million records.


texts:
Changing the “ArcGIS Server Data Extraction tool” default directory

I have successfully published the data extraction tool and its working properly, but I am wondering if I can configure the tool in order to change the output folder or not.

If I use the tool that the at the level of the desktop I can edit and choose the default folder

The published tool the output folder cannot be edited!

How Can I configure the published tool so that the user can control and change the output folder location.


texts:
I have a raster area which i could also convert into a vector area. I need all cells which lie in direct vicinity of the edge(s).

My first try was to convert the vector area to lines and then the resulting border line back to raster. The problem with this approach is that i won't get all border cells- but instead just some and spaces in between (see picture).

Any better idea how to solve my issue?


texts: Problem in storing splitted strings in attribute table

texts:
I want to split strings of column 'Text' and save them in column 'Text_edit'. So far the split and storing the new values in a new column works. However as can be seen in the attached figure, only the last value seems to be stored in every feature of the second column.
I aparently have a problem in my code for storing the strings in the new column.
Below my code:
 #example string to be splitted: {\fArial|b0|i0|c0|p34;\C6;A1 (ab GOK)}
 #split output => A1 (ab GOK)

        b.dataProvider().addAttributes([ QgsField("Text_edit", QVariant.String) ])
        b.startEditing()

        for f in b.getFeatures():
             split1 = feature['Text'].split(";")
             split2 = split1[-1]
             split3 = split2.split("}")
             split_final = split3[0]
             newval = split_final
             f.setAttribute(f.fieldNameIndex('Text_edit'), newval)
             b.updateFeature(f)

        b.updateFields()

        iface.vectorLayerTools().stopEditing(b,False)


texts: How to overlay WFS output in GML format on a googlemap api v3 using GeoXML3

texts:
I would like to overlay a polygon layer on Google map API V3  from a WFS output file  in GML2 format from a remote WFS server. No other output format is available in the WFS server.

Is it possible to do this using GeoXML3 ?

It would be great if some can give me a demo link  which using a GML2.


texts: Trying to activate GeoRover LocusTracker toolbar; another way to import to shapefile?

texts:
In ArcGIS I'm trying to activate the GeoRover LocusTracker toolbar but it won't come off its grayed out status. Might anyone know how to make it active?

If anyone knows of another way to import an Excel file or convert an Excel file into a shapefile (with just fields of FID numbers, province names, and MGRS).


texts: How to make a pgr_dijkstra wrapper with reverse cost?

texts:
I'm trying to prepare a dijkstra wrapper for Geoserver. I got a code from my professor and have to modify it so that it uses the dijkstra algorithm with the reverse cost argument. Can anyone please tell me what I'm doing wrong?

I get this error: Support for id,source,target columns only of type: integer. Support for Cost: double precision.
Here's the code:

CREATE OR REPLACE FUNCTION pgr_fromatob(
 IN sche character varying,
 IN tbl character varying,
 IN x1 double precision,
 IN y1 double precision,
 IN x2 double precision,
 IN y2 double precision,
 OUT seq integer,
 OUT gid integer,
 --OUT name text,
 OUT heading double precision,
 OUT cost double precision, --nie wiem jak podac cost i reverse_cost do record
 OUT cost_reverse double precision,
 OUT geom geometry)
 RETURNS SETOF record AS
$BODY$
DECLARE
 sql text;
 rec record;
 source integer;
 target integer;
 point integer;

BEGIN
-- Find nearest node
EXECUTE 'SELECT id::integer FROM osm_rowerowe_500_points_y
ORDER BY geom <-> ST_GeometryFromText(''POINT('
|| x1 || ' ' || y1 || ')'',2180) LIMIT 1' INTO rec;
source := rec.id;
EXECUTE 'SELECT id::integer FROM osm_rowerowe_500_points_y
ORDER BY geom <-> ST_GeometryFromText(''POINT('
|| x2 || ' ' || y2 || ')'',2180) LIMIT 1' INTO rec;
target := rec.id;

-- Shortest path query (TODO: limit extent by BBOX)
 seq := 0;
 sql := 'SELECT id, geom, czas, czas_reverse, source, target,
ST_Reverse(geom) AS flip_geom FROM ' ||
 'pgr_dijkstra(''SELECT id as id, source::int, target::int, '
 || 'czas::float AS cost, czas_reverse::float as reverse_cost FROM '
|| quote_ident(sche) ||'.'
|| quote_ident(tbl) || ''', '
 || source || ', ' || target
 || ' , false, true), '|| quote_ident(sche)
||'.'|| quote_ident(tbl) || ' WHERE id2 = id ORDER BY seq';

-- Remember start point
 point := source;
 FOR rec IN EXECUTE sql
 LOOP

-- Flip geometry (if required)
IF ( point != rec.source ) THEN
rec.geom := rec.flip_geom;
point := rec.source;
ELSE
point := rec.target;
END IF;

-- Calculate heading (simplified)
EXECUTE 'SELECT degrees( ST_Azimuth(
ST_StartPoint(''' || rec.geom::text || '''),
ST_EndPoint(''' || rec.geom::text || ''') ) )'
INTO heading;
-- Return record
 seq := seq + 1;
 gid := rec.id;
 --name := rec.typ_d;
 cost := rec.czas;
 cost_reverse := rec.czas_reverse;
 geom := rec.geom;
 RETURN NEXT;
 END LOOP;
 RETURN;
END;
$BODY$
LANGUAGE plpgsql VOLATILE


texts: Removing water from OpenStreetMap polygons?

texts:
I have a full OSM import using non-customized osm2pgsql (slim).  I need to remove all coastal water from the polygons requested by the user on the fly.

I have imported all water from osm water polygons shapes using shp2pgsql -I -s 3857 -g way water_polygons.shp water_polygons | psql -d gis, which creates a water_polygons table with 41,000+ rows.

It seems that ST_Difference is the PostGIS function to use, but I'm not sure of the best/most efficient approach.

Should I merge all the water rows into one row to speed up
ST_Difference, or subtract each water polygon one by one?
I combine several polygons into one. Should I subtract water before combining the polygon, or after?

This is my current code to extract GeoJSON polygons without water subtraction. The polygons are grouped by their Wikidata ID.

SELECT id,
 ST_AsGeoJSON(ST_Transform(ST_ForceRHR(way), 4326)) as data FROM
(
  SELECT id, ST_Multi(ST_Union(way)) AS way
  FROM (
    SELECT tags->'wikidata' AS id, (ST_Dump(way)).geom AS way
    FROM $1~
    WHERE tags ? 'wikidata' AND tags->'wikidata' IN ($2:csv)
    ) tbl1
  GROUP BY id
) tbl2


texts: MSSQL Spatial vector features not showing up in QGIS 2.18 up

texts:
I'm importing a MSSQL 2008 Spatial layer in QGIS. In version 2.14 it works OK, I can see the vector features and interact with the associated attribute tables. But in versions 2.18 and 2.99, though the connection is OK (I can see the attribute tables), the features are not showing up in the map canvas. I can even modify the data in the tables but I can't graphically select the features for editing.

I'm using a Windows 10 PC. The 3 versions of QGIS are runnig in the same machine, at the same time, with the same configuration. I've tried with the standalone installer and with the OSGEO4W installer (32 bit).

¿Anything changed from version 2.14 to newer versions that I can't get the features drawn on the map?

ADDED 2017-10-23

In QGIS 2.18.13, in the Layer Properties dialog for the CAMPOS_geo layer, under Metadata, there's a Properties section with this data:

General
    Storage type of this layer
    MSSQL spatial database

    Description of this provider
    MSSQL spatial data provider

    Source for this layer
    dbname='GIS' host=SERVER2 user='sa' srid=4326 type=Point table="dbo"."CAMPOS_geo" (CASCO_point) sql=

    Geometry type of the features in this layer
    Point (WKB type: "Point")

    The number of features in this layer
    611

    Capabilities of this layer
    Add Features, Delete Features, Change Attribute Values, Add Attributes, Delete Attributes, Create Spatial Index, Create Attribute Indexes, Fast Access to Features at ID, Change Geometries

Extents
    In layer spatial reference system units
    xMin,yMin -65.4273,-31.905 : xMax,yMax 0,0

    Layer Spatial Reference System
    +proj=longlat +datum=WGS84 +no_defs

And in this picture is the configuration of the CAMPOS_geo table in SQL Server.

QGIS 2.14.19: Vector Features and associated data shown OK.

QGIS 2.18.13: Data shown but Vector Features not shown.

QGIS 2.99.0: Data shown but Vector Features not shown.

 MSSQL Spatial table setup


texts: CreateEnterpriseGeodatabase_management fails on PostgreSQL claiming ST_GEOMETRY not present

texts:
I have a PostgreSQL server (9.3.10) already setup for ArcGIS enterprise GDBs with one GDB already in it (created in the past, with ArcGIS 10.3 I believe). I am using ArcGIS 10.5.1 on my client to create another GDB on it and it fails with the below message.

  You must copy the latest ST_GEOMETRY library to the PostgreSQL
  software location. Refer to the ArcGIS help topics for more details.
  Connected RDBMS instance is not setup for Esri spatial type
  configuration. Failed to execute (CreateEnterpriseGeodatabase).

Is this a version mismatch issue or bug? What is the work around besides upgrading the PostgreSQL instance (if that is a solution at all).


texts: Applying a Single Gradient Fill Color to a multipart feature

texts:
So I have an indicatrix covering a world map for the purpose of illustrating distortions for different projections. I want the circles of the indicatrix to share the same gradient fill color, as opposed to each having their own individual gradient fill. To clarify, here is a picture of the indicatrix:

Notice that each circle has its own individual gradient fill. What I want to do is have QGIS treat them all as part of the same shape, so that a single gradient fill is spread out across them all. To visualize this, take the gradient fill for this circle:

...and then use the indicatrix to "clip" the circle without actually changing the fill. (I literally tried this already fyi. Just creates an visually identical shapefile to the indicatrix.)

I figured all I had to do was dissolve my indicatrix shapefile to a single feature and it would replicate this result, but apparently the issue is more complex than that. Can anyone enlighten me?


texts: Creating line (closest vertex to line) using ArcGIS Desktop?

texts:
I'm using ArcInfo 10 SP3.

I am working in reorganising our utility data.  Two years ago, we have started collecting  Private Service water line.  We still have a lot of them to extract from old Record Drawings.

I was wondering if there was a way to create line that would join our building footprints to the WaterMain line?

I would like to use the the building vertex closest to the Water main as a starting point.


texts: Researching into the nuts and bolts of displaying ESRI Grids in QGIS

texts:
I am trying to load an Arc/Info Binary Grid into QGIS 2.2.0 Valmiera installed in OS Windows 8.0. As far as I know, I only need to load the "w001001.adf" file contained in the "altitud" folder while keeping all other files there (Fig. 1), to see the raster and be able to work with it. There is no files names "altitud" to be loaded (as in ArcMap), and I have to load the "w001001.adf", which appears under the very same name in the layers toolbar (Fig. 2).

Conversely, when using ArcMap, we start adding the same raster data from the very same folder, and only one file with the name of the layer (altitud) appears (Fig. 3). All other files (dblbnd.adf, hdr.adf, prj.adf, sta.adf, w001001.adf, and w001001x.adf) are not displayed. One only needs to double click on the "altitud" file to load it, and this layer (and others) is displayed under the same name in ArcMap (Fig. 4).

At this point my questions are:

1) Why does QGIS behave somewhat differently from ArcMap in this respect?

2) Is there a way to load and see raster files in QGIS in the same was as in ArcMap, getting to see the actual layer names in the layers toolbar?


texts: Creating web map using MapGuide Maestro?

texts:
I have been given a task by my lecturer to create a web map using mapguide maestro or Autodesk Infrastructural Studio. The application that I used to finish my mapping are mapguide maestro, and the data is in shapefile. My mapguide open source are still installed and running.

Is there any site that has a tutorial to do it? I am trying to make a website that looks like this.


texts: Changing pixel size from raster by value

texts:
How to change pixel size from raster by value using PyQGIS or Python GDAL?

I tried to using SAGA GIS Tool Resampling:

ext2 = rasterLyr2.extent()
xmin = ext2.xMinimum()
xmax = ext2.xMaximum()
ymin = ext2.yMinimum()
ymax = ext2.yMaximum()
extends= "%f,%f,%f,%f" % (xmin, xmax, ymin, ymax)
processing.runalg('saga:resampling', rasterLyr2, True, 0, 0, extends, 10.0,None)

But that way not work for me because in the some CRS have different extends order any idea how to change pixel size from raster using PyQGIS?


texts: Arrows with direction and proportional to attribute

texts:
my data looks like this:

> main
   end.station start.staion Freq end.latitude end.longitude        end.station
1          359          519  929     40.75510     -73.97499 E 47 St &\nPark Av
2          477          465 5032     40.75641     -73.99003   W 41 St &\n8 Ave
3          484          519 1246     40.75500     -73.98014   W 44 St &\n5 Ave
4          484          318 2654     40.75500     -73.98014   W 44 St &\n5 Ave
5          492          267 1828     40.75020     -73.99093   W 33 St &\n7 Ave
6          492          498  957     40.75020     -73.99093   W 33 St &\n7 Ave
7          492          362 1405     40.75020     -73.99093   W 33 St &\n7 Ave
8          493          477 1582     40.75680     -73.98291   W 45 St &\n6 Ave
9          493          529  728     40.75680     -73.98291   W 45 St &\n6 Ave
10         529         2021 1748     40.75757     -73.99099   W 42 St &\n8 Ave
   start.latittude start.longitude              start.station
1         40.75188       -73.97770         Pershing\nSquare N
2         40.75514       -73.98658        Broadway &\nW 41 St
3         40.75188       -73.97770         Pershing\nSquare N
4         40.75320       -73.97799 E 43 St &\nVanderbilt\nAve
5         40.75098       -73.98765        Broadway &\nW 36 St
6         40.74855       -73.98808        Broadway &\nW 32 St
7         40.75173       -73.98754        Broadway &\nW 37 St
8         40.75641       -73.99003           W 41 St &\n8 Ave
9         40.75757       -73.99099           W 42 St &\n8 Ave
10        40.75929       -73.98860           W 45 St &\n8 Ave

I've managed to connect the end.station and start.station points in QGIS using the MMQGIS Hub lines tool.   However, I am unable transform these lines into arrows that a) point in the direction of the path, and b) have a width that is proportional to FREQ.

I have tried with the Arrows plugin, but to no avail.  All I get are very long lines that extend infinitely.  How does one solve this?


texts: Creating random points using density-QGIS

texts:
I'm trying to generate random points at a specific density within a polygon layer.  The polygon layer has 1850 features and I'm looking to have random points within each feature at 1/meter. Going through Vector -> Research tool -> Random points -> "Use this density of points" = 1 (assuming the default is based on unit of the project) -> Save As: .shp -> click OK.  The OK button darkens, then lightens to normal without producing any result.

I've also tried going through the Processing Toolbox -> QGIS geoalgo -> Vect Creation Tools -> Random points inside polygons (fixed), set density = 1, and get error:

"Algorithm Random points inside polygons (fixed) starting...

float division by zero See log for more details"

Any ideas? I think it has something to do with the density parameter, but I haven't found what the unit used for the density is.  Insight welcomed.

QGIS 2.14.3-Essen


texts: Sorting text column by numbers with letters

texts:
I have a text column. It's made up of the numbers 1, 2, 3 up to 1650, which are currently not in increasing order. However, it also includes text such as 520a, 639a, 639b, 639c. Sorting on increasing value would produce 1, 10, 100, 1000, etc, rather than 1, 2, 3, ...
I know how to convert from text to integer, which would work if my cells did not include alpha characters.
How can I sort to end up with 1, 2, 3, ..., 10, ..., 100, ..., 519, 520a, 521, ..., 638, 639a, 639b, 639c, 640, ..., 1650?


texts: I've been given a raster with flood elevation limits along a river, it's not a constant flood level so elevations range by ~35m.  I need to do an analysis of the impacts if these elevations were 2m higher - I know I can raise the elevations in the raster calculator, but how can I project where the new limits would fall horizontally on the underlaying surface raster?


texts:
I have accomplished this task in the field calculator (thanks to another post here, but I now realize that I need to do it using arcpy).

I have code that perfectly marks the second, third etc. duplicate(s) but I need it to mark all values that have a duplicate.  Basically, I can't figure out how to mark the FIRST instance as well as the rest.  Here's the code:

import arcpy

inFeatures = "thephaseIIIlayer"
fieldName = "PhaseIII_ID"

desc = arcpy.Describe(inFeatures)
fields = desc.fields
rows = arcpy.UpdateCursor("thephaseIIIlayer", "", "", "PhaseIII_ID", "PhaseIII_ID A")

    i = 0
    for row in rows:
        if i == 0:
        value = row.PhaseIII_ID
    i += 1
    if i > 1
        row.SplitMarker = i
    if row.PhaseIII_ID != value:
        value = row.PhaseIII_ID
        i = 1
    row.SplitMarker = i -1
    rows.updateRow(row)

del row
del rows


texts: Querying ArcGIS MapServer layer in QGIS?

texts: I can add an ArcGIS map layer to my QGIS 2.18 project, or use the existing project provided by GeoNB (http://geonb.snb.ca/documents/misc/geonb_qgs.zip).  Using the identify tool, I can see features of different entities in each layer. For example, I can select a parcel from the parcels layer and it returns, among other things, a PID and an area.  I would like to be able to search within these layers to find the geometry I am looking for.  For example, I would like to be able to, knowing a PID, locate a geometry in my viewport.


texts:
How can I use Leaflet and OpenStreetMap to display 3D buildings for free? I find same scripts but needs API key or access token. Who already has a script to show 3D map Leaflet using open source?

i tried with osmBuilding but the codes which i find doesn't works good, have scripte .hmtl work well?


texts: Generating convincing details to low resolution DEM

texts:
For non-serious fields, game-terrain editing to be specific, I've got to generate much more detailed DEM (5-meter at least) from 1-arc-sec SRTM, which has the best resolution I can find yet.

I wonder if there's any interpolation / noise-generation / etc. algorithm that can re-sample lower-resolution DEM and create details accordingly without harming the credibility of the terrain to human eyes.


texts:
I'm assisting a retired user in combining a large number of raster tiles downloaded from Historic Digimap http://edina.ac.uk/digimap/description/historic_overview.shtml into one project in qgis.

I'd like if possible to combine these into one layer, - I'm aware that these can likely be stitched together using an external tool, but am unsure exactly how to go about doing so (or whether this is in fact possible in qgis itself) in a way that the resulting layer is still usable on a standard desktop pc (~4GB ram) or whether this is in fact not possible.

Many thanks in advance


texts: Loading Scene from ID into ArcGIS JavaScript API 4: "SyntaxError: Unexpected token ILLEGAL"

texts:
I have created a city engine web scene,and I am trying to display the 3d-visualisation view on sceneview by using the portal Item Id of my webscene. But I am getting

  SyntaxError: Unexpected token ILLEGAL

How to solve this?

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
  <title>Load a basic WebScene - 4.0beta3</title>

  <link rel="stylesheet" href="https://js.arcgis.com/4.0beta3/esri/css/main.css">


texts: ArcMap exporting nonsensical data to shapefile

texts:
I am using ArcMap 10.5.

I am trying to make a copy of a feature class to make some edits to it. The feature class I am trying to make a copy of is a File Geodatabase Feature Class, and I am trying to export it as a shapefile by clicking 'Data -> Export Data...' from the table of contents. This creates an output file without any errors, but when I add the new exported data back to the map and open the attribute table, I just get garbled mess.

A portion of the original attribute table:

The same fields of the copied attribute table:

No idea why this is happening as I'm not getting an error message and it seems a fairly routine operation that I'm pretty sure I have carried out before... I note as well that not all of the rows have been exported.


texts: Converting NZTM (New Zealand Transverse Mercator) to Lat/Long

texts:
This link is pretty close to what I am aftering (Converting NZMG (or NZTM) to latitude / longitude for use with R map library).

Basically I am looking for a formula that can convert NZTM coordinates to lat/long (to be run on C# platform).  I need this implemented in a C# code that can take in a pair of NZTM coordinates and return lat/long result. In other words:

public double[] (double nztmLat, double nztmLong){
//Implementation details
//return lat/long in a double array.
}


texts: How to extract X Y coordinates and RGB from pixels using python in QGIS?

texts:
I'm trying to white a Python plugin in QGIS (2.2) to get, for each pixel of a raster image, its X Y coordinates, as well as, its RGB (3 bands) values.

At first I opened the raster file using: rlayer = QgsRasterLayer(rfileName, rbaseName)

Now I don't know how to get, for example, for pixel (1,1) its coordinates (X,Y) and its RGB color values.

Anybody could help me?

I know that I'll need to implement a DO WHILE, but I don't know the commands to extract this information of each pixel.

Thanks in advance,
Mateus


texts: Error accessing spatial database connection with ArcGIS 10

texts:
I set up a new SDE database with ArcGIS 10 and everything appeared to go correctly.  The service is running and I'm able to view the database in SQL Server.

However, when I try to access the spatial database connection in ArcCatalog 10, I get the following error:

  Failed to connection to the specified
  server. This release of the
  GeoDatabase is either invalid or out
  of date. DBMS table not found
  [Microsoft SQL Server Native Client
  10.0: Invalid object name 'sde_mymap.dbo.GDB_Release'.]
  [sde_mymap.dbo.GDB_Release]

Both ArcSDE and ArcCatalog are version 10.

I checked the sde_setup.log file and I can find no record of that table being created, or any other errors.

Please advise.

Thanks,
Seth


texts:
Is there a tool (or plugin) to extend a segment in Qgis (3)?

I would use it to find objects on a map which are visible on photographs taken from known locations but too far away from those locations to simply search on the map by following a rough direction. So if I draw a segment from the viewpoint to a close object (easier to find on the map) which is in the same direction as the far object and extend that segment, I can then navigate the map knowing the far object is somewhere along that line.


texts: Sharing large raster image with others on webservice/stand alone webpage

texts:
I would like to share part of the NOAA/NGDC VIIRS Nighttime imagery (found at https://ngdc.noaa.gov/eog/viirs/download_monthly.html) with other people within my organization that do not use a spatial analysis software package. In my case this is QGIS. What I would like to do is upload the map I have created to some cloud service (obvious candidate here would be QGIS cloud), or share it with other people while the data is embedded in some sort of html file. Let's say I would like to share the Nighttime imagery of one country only, which we'll take as Ivory Coast for simplicity's sake. Here are the steps I have undertaken so far:

Downloaded the latest nighttime light image for November 2016. To be able to crop the image to the right size, I have also downloaded the natural earth dataset and used the ne_10m_admin_0_countries vector layer.
Once both data sets have been loaded into QGIS, I have cropped the nighttime imagery with the Raster>Extraction>Clipper tool. Hereby I use the Ivory Coast shape as the mask layer. To reduce the size of the 'clipped' nighttime imagery I also use suggestions found here: Clipping raster-image increases file size. (I fear that my troubles start here, since the 'clipped' image still contains the entire extent of the original nighttime image, even though I only want to have the Ivory Coast part. This manifests itself as everything being black outside of the cropped area up to the extent of the original image. I have looked for possibilities of 'cutting'/
'cropping' (as you would do in any image manipulation software) all of the other data off, but couldn't find anything)
Export the obtained results with some kind of tool. So far I have tried:

QGIS2Web: When I try this, it takes a long time before I get to the dialogue window where I can select relevant options. But as soon as I select the 'leaflet' option (the other option just shows a white screen), QGIS locks up and I am forced to restart the program. Retrying this does not do anything.
QGIS Cloud: When I try to upload the cropped raster layer, the uploading is extremely slow since it never comes close to the number of tiles it needs to upload (or maybe I am mis-assuming the size of the data that needs to be uploaded). In any case QGIS hangs often with this option and crashes after a while.

The other option to share these maps would be to just print them out, but the map maker also crashes if I want to export to an image/pdf file.

Could this be an issue with my PC (i5, 16GB RAM, Windows7), or am I trying to push too much data through QGIS and it's extensions?


texts:
Do you know any good tutorials or any resources for learning from the beginning GML?

I am asking  because (later) i would like to learn how to make a application schemas in context of INSPIRE.


texts: GeoServer: How do I get rid of error Unable to parse number For input string: "x"?

texts:
I am using GeoServer version 2.17.2 at my Windows 10 (x64) machine. In the PostGIS store, I have an SLD symbology coupled layer. Using GeoWebCache, I have created tiles of the same, and I could see several PNGs in Geoserver data directory (Windows machine root directory):

At the moment, I am trying to display (without any mapping client, OL3, Leaflet, etc.) these tiles (via GWC TMS service) in my Google Chrome browser using:
http://my_server:my_port/geoserver/gwc/service/tms/1.0.0/my_workspace%3Amy_layer@EPSG%3A900913@png/{z}/{x}/{y}.png

But, I keep getting the error:
400: Unable to parse number For input string: "{x}" from /service/tms/1.0.0/my_workspace:my_layer@EPSG:900913@png/{z}/{x}/{y}.png

I tried to google but couldn't get any hint to resolve the above error. Can someone help me to understand what am I doing wrong here?


texts: Making transparent layers with Google Earth layer

texts:
I managed to get Google Earth layer, but when I add my GPS gpx files, all are hidden under the Google map...
How to put them onto Google Earth layer?


texts: Using UpdateCursor to update a field based on other field values

texts:
I've been trying to create a python script that will update certain fields in the arcgis attribute table. I have attempted to do this by using the update cursor to update a field based on another field value found in the same row.

Here is the code I used:

import arcpy

fc = "C:\Users\maureen\Documents\ArcGIS\EDRN\EDRN_LINK.shp"
fields = ["FID", "RIVERNAME"]

with arcpy.da.UpdateCursor(fc, fields) as cursor:
    for row in cursor:
        if  row[0] = 4401:
            row[1] = "Aire"
            cursor.updateRow(row)

Unfortunately, when I run this in ArcMap it returns the following syntax error:

Parsing error SyntaxError: invalid syntax (line 9)

Does anyone have any idea where I'm going wrong? I've researched this script a fair bit and can't see any obvious errors.


texts: React-Leaflet: Marker has wrong position after zoom out

texts:
When I show a marker in the center of the map it is OK.
But when I zoom out, the marker is like staying at the same place. Finally its position is very wrong.

Why is it not moving together with the map when zooming out?

I user react-leaflet:

const customMarker = new L.icon({
    iconUrl: '/images/icons/location-pointer.svg',
    iconSize: new L.Point(35, 46),
});

...

renderMarker(position) {
    if (position) {
        return (
            <Marker
                position={[position.lat, position.lng]}
                icon={customMarker}
            >
            </Marker>
        );
    }
}

...

render() {
        return (
            <Map
                ref={this.mapRef}
                center={positionMap}
                zoom={position.zoom}
                onClick={this.handleClick}
                onMoveend={this.handleMoveEnd}
                dragging={true}
                scrollWheelZoom={false}
            >
                <TileLayer
                    attribution="&amp;copy OpenStreetMap contributors"
                    url="https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}"
                />
                    { this.renderMarker(this.props.locationData.locations[0]) }

            </Map>
        );
    }
}


texts: Update joined table fields when a MySQL database is updated

texts: Using Qgis I joined feature tables with tables from a MySQL database. I imported the database's tables as new vector layers. Is there any way to update the joined table's fields each time the MySQL database is updated?


texts: Under what situation do we use Rasterisation or Vectorisation?

texts:
As mention in the title, Rasterisation is the conversion of vector maps to raster grid. While Vectorisation is the conversion of raster grids to vector maps but under what situation would we use these?

Raster maps are better at continuous data,
Vector maps are better at discreet data.

What is the point Rasterisation or Vectorisation, for example if raster map are better at continuous data but you are trying to force it into a vector map.


texts: How to drag a layer from ArcMap Add-in into its MXD?

texts: I have a simple ArcMap Add-in that displays all featureclass names in a database inside of a ListBox.  If I select a name or multiple names in the listbox, how can I drag and drop it in ArcMap?  Just to give a quick summary, the purpose of this add-in is because the database has encrypted names and I have to show the "friendly names" in the ListBox control so I will use a special refrence table in SQL to populate the listbox with the friendly names.


texts: Detecting if point is on left or right side of line in PostGIS?

texts:
I have a linestring table and a point table in postgis.

I know the closest line to any given point.
What I need to know is on which "side" of that line is the point.
I guess I have to do that by creating a perpendicular line from given point to the line (closest point on the line) and then compare the coordinates, but I don't know exactly how to do that, and if it's the proper way, since line changes it's direction.

I've made a picture to illustrate my task.

The line itself is black, its direction is shown with green arrows.
I need to add a "side" column to the point table, so that red points should have value "right" and blue points should have value "left".

Can someone give an SQL code example of calculating a "side" value of a point?


texts: Control colour point on ramped gradient in QGIS symbology based on value

texts:
Is there a non manual way to control the colour ramp in symbology based on value.

In my example below I have values that can be positive or negative, i want everything positive to be one colour and negative the other with values around zero as white.

How does QGIS decide where to place the break in colour? Double clicking on the ramp brings me this set of options but the figures at bottom don't correspond to the figures in the source data.

Answer given by Eric below, screen grab is that applied to my example for those who can't translate the German. This still leaves you with some manual edits to make, for most automated solution see my solution


texts: Understanding extent used by ArcGIS Create Fishnet?

texts: I am trying to create a 20 X 20m Fishnet for the island of Wallis. But I am having an issue with the output. As you can see in the image below the grid is not created and the extent is way out of the island polygon shapefile. I have done it for my own country without any issues but every time I do it for this island the fishnet does not follow the proper extent that I set out for the island.


texts: How could you calculate the distance between multiple associated points, specifically a list of structures serving specific addresses within CartoDB?

texts:
I have loaded two point datasets into CartoDB. The first is a list of structures serving specific addresses, each structure has a unique code associated to it as well as the location of the structure itself. The other point dataset is a list of the addresses and their locations, each address having the code of the structure that serves it.

I know how to calculate the distance of each structure from each address, however this seems quite inefficient when all I require is the distance between the structure serving the address in question (the other structures I'm not worried about for any given address).

Note I only require the great circle (haversine) distance. Also, if the same query can determine which of the addresses lie X distance (in km) away from the structure serving it that would be great.

Any ideas on how you'd go about this within CartoDB/PostgreSQL?


texts: Why is the data missing in my new shapefile?

texts:
I've made a new shapefile from a part of an attribute table. Marked all the rows that I needed and clicked on "save selection as", added the new layer in my project. All the data is missing, there is only the headlines of the columns.

I save the new shapefile like this:

Format: ESRI shapefile
Encoding: System
Layer CRS
And the same projection as all the other layers, SWEREF99 12 00.

I'm used to working in ArcMap, and I'm now working in QGIS. Is there a difference in the programs of how I do this kind of thing?


texts: Openlayers text symbology is not being shown

texts:
I have a Vector Layer in a .json file:

{
    "type" : "FeatureCollection",
    "features": [
        { "type": "Feature", "properties": { "Station_Id":107, "NAME":"ARDEBIL(F)","ELEVATION":"1314.3"}, "geometry": { "type": "Point", "coordinates": [ 48.24 , 38.19 ] } }
    ]
}

For symbology I use an SLD, and the SLD File is following:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
  xmlns="http://www.opengis.net/sld"
  xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>me</Name>
    <UserStyle>
      <Title>GeoServer SLD Cook Book: Rotated square</Title>
      <FeatureTypeStyle>
        <Rule>
          <Name>salam Bar Mahdi 2</Name>
          <Title>Title 2</Title>
          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>circle</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#00FFFF</CssParameter>
                  <CssParameter name="fill-opacity">1.0</CssParameter>
                </Fill>
                <Stroke>
                 <CssParameter name="stroke">#000000</CssParameter>
                 <CssParameter name="stroke-width">2</CssParameter>
                 <CssParameter name="stroke-opacity">1.0</CssParameter>
               </Stroke>
             </Mark>
             <Size>6</Size>
             <Rotation>0</Rotation>
           </Graphic>
         </PointSymbolizer>
         <TextSymbolizer>
           <Label>
             <ogc:PropertyName>NAME</ogc:PropertyName>
           </Label>
           <Font>
             <CssParameter name="font-family">Arial</CssParameter>
             <CssParameter name="font-size">120</CssParameter>
             <CssParameter name="font-style">normal</CssParameter>
             <CssParameter name="font-weight">bold</CssParameter>
           </Font>
           <LabelPlacement>
             <PointPlacement>
               <AnchorPoint>
                 <AnchorPointX>0.5</AnchorPointX>
                 <AnchorPointY>0.0</AnchorPointY>
               </AnchorPoint>
             </PointPlacement>
           </LabelPlacement>
           <Fill>
             <CssParameter name="fill">#000000</CssParameter>
           </Fill>
         </TextSymbolizer>
       </Rule>
     </FeatureTypeStyle>
   </UserStyle>
 </NamedLayer>
</StyledLayerDescriptor>

I can see only a circle symbology and I can't see the text symbology.


texts: Join was successful however resulting table has no rows, create table with now rows?

texts:
I am trying to spatially join a layer of points to a polygon layer using ArcGIS Desktop 10.4. Both layers have the same projection system.

I am getting an error message:

  The join was successful however the resulting table has no rows,
  create table with now rows?

Then an empty layer is created, that has the correct column headings, but no rows.

Anyone know what the problem could be?


texts: How to calculate a curve distance between two xyz points?

texts:
I have a curve surface (example is attached).
Each node has its own coordinate in geographic system (longitude,latitude and depth). I want to calculate a distance between two nodes (from the one on the upper left to another on the lower right). However, I want it a distance that follows curve of the surface (not a shortest distance).

Any idea how to do that in python?


texts: How to display the label with "Identify Features" in QGIS?

texts:
I have been struggling for one hour with this simple problem (QGIS 1.8.0-Lisboa, Windows 7) :

I have a raster layer. I would like to display the label and not only the band value when I use the Identify Features button. Is it possible?


texts:
So, I'm trying to draw curved lines in my application which is based on a Leaflet map.

As far as I know Leaflet does not support drawing curved lines at the moment, so my question is really more about how to best get some curved line funtionality in there.

The problem is made more difficult because my code will have to be IE-8 compatible, which means that pure SVG is no solution. (Leaflet actually does a fall-back to VML if it detects that no SVG support is present...)

So one possibility would be to myself write some leaflet extension code based on SVG with its own fallback to VML. This would be a hell of a lot of work. :/

Does anyone have a better proposition?


texts: Does Oracle Spatial calculate geodetic distance?

texts:
I've got to migrate a PostGIS database to Oracle Spatial, and am stuck on a problem:
In a view, I have a column that amounts to ST_Distance(geography(geom_1), geography(geom_2)), effectively calculating geodetic distance in meters.
In Oracle, calculating the distance in WGS-84 gets me the result in degrees. I've also tried to project the geometries to World Mercator and calculate distance, but when data gets far from the equator, the results become more and more different.
How can I replicate this code in Oracle Spatial?


texts: Dynamic filtering based on selected Layer in ArcGIS Desktop?

texts:
I am not sure this is even possible, I have about 38 geographical polygon layers showing areas of the UK, I also have a query (point layer) that shows all of the offices around the uk. What I am attempting to do is if I select layer 1 that only the offices that intersect with displayed polygons are show. Then if I display an additional layer or change the displayed layer the shown points change accordingly.
Image one shows offices outside and inside the layer polygons,

Image two I have switch the polygon layers but I want only the intersecting offices to show.
I could do this with having 38 separate attribute queries but that seems an excessive amount of work.

I am using ArcGIS Desktop 10.1.


texts: Maintain 8-bit pixel depth in Raster Calculator

texts:
Would anyone know how to maintain 8-bit pixels when computing raster calculations in ArcMap?

This is part of a model so I need to keep the image with the original 8 bit, do the calculations, convert it into 32bit and do some other calculations and then output it.

The problem is to keep it at 8 bit for the first set of calculations.


texts: Using a local WMTS server with ArcGIS Earth on iPad

texts:
I would like to use an open map tile server to server tiles to the ArcGIS Earth application running on my iPad. I would like to do this as the iPad is not always connected to the Internet, but will always be on the same network as the tile server.
Currently the tile server is running on a Synology NAS using the klokantech-tilserver-gl docker package. The ArcGIS Earth iOS version is 1.3.0
I tired to add a data source by adding the url for the tile server, an example is

http://192.168.1.32:8080/styles/klokantech-basic/wmts.xml

Also tried

http://192.168.1.32:8080/styles/klokantech-basic/wmts
http://192.168.1.32:8080/styles/klokantech-basic/

Always get the same error messgae "Please input correct URL".


texts: Check if Point belongs to LineString in Geotools

texts:
Is there any way to check with Geotools if a point is placed on a LineString. I am trying to retrieve street names of a route obtained from a routing algorithm. The route is displayed in gps coordinates and i was thinking to look for the name of the street by checking if the point is placed on the LineString.

If you have other suggestions, i would be grateful.


texts:
I have a point feature that I would like to buffer to select segments that are within 15 meters of the point. The code I have is working but it's influenced by the map scale. Zooming in or out changes the buffer size. This is what I am using.

pntGeom = feature.geometry()
pntBuf = pntGeom.buffer(canvas.mapUnitsPerPixel()*15, 0)

is there a way to consistently buffer the point regardless of zoom level?


texts: Relation between scale factor and standard meridians/parallels

texts:
Let's use transverse Mercator with a custom zone, with central meridian 30, the zone is 3 degrees wide at the equator, with a scale factor k=1, it means that the cylinder touching the spheroid along the central meridian.

However, with scale factor k=0.9996, the cylinder is secant of the spheroid at two meridians, one at right of Central Meridian and the other at left. How to find these meridians for a given scale factor k at equator?


texts: If Else Statement Existence of a File With User Input

texts:
I am trying to design a block of code that checks if a user inputted directory exists, prints if it does, and if it does not.

 import os

 foldername = raw_input ("Please enter the name of your desired folder: ")

 if os.path.exists(foldername):
      print "Yes, your folder exists!"
 else:
      print "Please re-enter your folder directory as it does not exist."

The abov is what I have, but it consistently outputs that the directory does not exists, despite me knowing that it does.

Insights?


texts:
For my Google Earth Engine project, I am finding the mean NDVI for each year between 2000 and 2020 in protected areas a part of the World Database on Protected Areas (WDPA). As an example, I have included code written to find the yearly mean NDVI values in Salonga National Park, which is located in the Democratic Republic of the Congo (DRC). I am using the MODIS vegetation index produced on 16-day intervals.
The code finds the mean NDVI value within the park boundaries and then places it into a dictionary with the years as keys and mean NDVI as values. Is there an easier/more concise way to write this code? This example finds the mean NDVI for 2000 and 2001 so far.
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');
var drcongo = countries.filter(ee.Filter.eq('country_co', 'CG'));
var styleParams = {
  fillColor: 'b5ffb4',
  color: '868686',
  width: 2.5,
  opacity: 0.6,
};
//Map.addLayer(drcongo, styleParams, 'DRC');
Map.setCenter(23.5, -2.5, 6);

var salonga = ee.FeatureCollection('WCMC/WDPA/current/polygons')
.filter(ee.Filter.and(
    ee.Filter.bounds(drcongo.geometry()),
    ee.Filter.eq('ORIG_NAME', 'Parc national de la Salonga')
  ));
//Map.addLayer(dataset, {}, 'Parc national de la Salonga');

var modisNDVI = ee.ImageCollection('MODIS/006/MOD13Q1')
              .select('NDVI')
              .map(
                      function(img) {
                        var rescaled_NDVI = img.select('NDVI')
                                                .multiply(0.0001)
                                                .rename('NDVI_rescaled');
                        return img.addBands(rescaled_NDVI);

  }
);

var modisNDVI_2000 = modisNDVI.filter(ee.Filter.date('2000-01-01', '2000-12-31'));
var modisNDVI_2000_clip = modisNDVI_2000.mean().clip(salonga);
var mean_dict_2000 = modisNDVI_2000_clip.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: salonga.geometry(),
  scale: 500,
});

var modisNDVI_2001 = modisNDVI.filter(ee.Filter.date('2001-01-01', '2001-12-31'));
var modisNDVI_2001_clip = modisNDVI_2001.mean().clip(salonga);
var mean_dict_2001 = modisNDVI_2001_clip.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: salonga.geometry(),
  scale: 500,
});

var dictionary = ee.Dictionary({
  2000: mean_dict_2000.get('NDVI_rescaled'),
  2001: mean_dict_2001.get('NDVI_rescaled')
});
print('Parc national de la Salonga Annual NDVI', dictionary);


texts: With CQL_FILTER on GeoServer 2.11.x, how do I do an ORDER BY?

texts:
I need to sort my output because I want to animate my output in Leaflet. I realize that I could sort the resulting GeoJSON in my JavaScript code but that is less efficient.

When I Googled it I got to the GeoServer 2.12.x documentation that suggested I use an addition to my query SortBy=&.  So here is what my URL looks like:

http://localhost:8080/geoserver/wfs?service=WFS&version=2.11.2&request=GetFeature&typeNames=hurricane:hurricane_point&PROPERTYNAME=(point,strength,category)&CQL_FILTER=storm_id=%27al122017%27%20and%20is_prediction=false&sortBy=(occur_date)&outputFormat=JSON

So that didn't work, no error so I assume that GeoServer ignored the SortBy.

Any ideas?

EDIT:  OK, I'm an idiot!  It was working all along, I just didn't recognize it.  I was judging whether it sorted correctly by the display order of the markers in Leaflet.  As user30184 pointed out, without including the SortBy field(s) in the output it's hard to ascertain if the sort worked correctly.  I did some experimentation with adding +Desc to the sort field to get it to sort backwards and did some key comparisons and see that all's well.  Adding credit because user30184 helped with correct syntax, plus I learned a bit from the answer.