texts:
I am really new to QGIS and I've imported this shapefile of UK with weird lines (the ones not denoted with the red crosses).

What is the fastest way to remove these lines so that it leaves me with just the boundary of UK?


texts: Storing Central European data, which SRS/which data type?

texts:
I am creating a PostGIS 2.X data base to store some 500 or so sites as points in it, which are distributed in the northeastern part of Germany, northern Poland and parts of the southern Baltic coast. I will also store some lines and polygons from the same in it, which will be needed to calculate distances and areas. In general an error of about 10-20 m would be acceptable for my task.

I was thinking of storing the spatial data as geography in ETRS89 Lat/Lon (EPSG 4258), which should allow for best compatibility over the entire area, but cost me some functionality compared to geometry. Since I am new to PostGIS, I don't know how strongly this will affect me.

So here are my to questions:
1) Is there a considerable downside to using lat/lon-data and geography for this task?
2) Is there a projected SRS I could use for this area without getting to much of an error to my measurements?


texts: Editing the description in the dashboard view of a Carto map once the map has been published?

texts:
It appears that once you write a description in the dashboard view of a Carto map (where it says "Add description...") the text is un-editable. Is there a way to edit the text?

The map on the left shows before text is entered, the map on the right shows after text has been entered and cannot be edited.


texts: Exporting Landsat 7 Surface Reflectance through Google Earth Engine

texts:
I'm trying to export Surface Reflectance data from Landsat 7 through Google Earth Engine. However, the dimensions don't seem to match up: previously with NVDI data, I printed the image collection size divided by the number of snapshots, temporal resolution of 8 days, per year to check that it had 16 years worth of data.

imgcoll = ee.ImageCollection('LANDSAT/LE07/C01/T1_8DAY_NDVI') \
    .filterDate('2002-12-31','2018-12-31')

imgcoll.size().getInfo() / 46
16

Trying to print the image collection size of Surface Reflectance data yields a confusing result. The wikipedia page on Landsat 7 indicates a temporal resolution of 16 days, which should mean that when dividing by 23, the result should be 16 years (2003-2018) as well.

imgcoll = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR') \
    .filterDate('2002-12-31','2003-12-31')

imgcoll.size().getInfo() / 23
64430.782608695656

Why is that the case, and what is the proper way of exporting Landsat 7 Surface Reflectance data? In addition, has the Landsat 7 dataset from Google Earth Engine been corrected for the SLC failure? This was not a problem encountered with exporting MODIS data.


texts: How to replace corrupt representations stored in geodatabase using symbology from layer file?

texts:
I have a lyr file with symbology that I want to store in the database. Symbolology is already stored in the geodatabase but the symbology levels has been corrupted migrating from 9.3 to 10. Thus, I want to replace it with the symbology from the lyr file.

How do I do that?

Symbology is stored in the representations of my feature class.


texts: Weighted Overlay Tool does not display all classes

texts:
In ArcMap 10.1 I have reclassified several layers into 4 classes according to suitability. I have then run the Weighted Overlay tool, but the output does only have classes 1,2 and 3. Does anyone have an idea why that might be? What happened to my class 4?

I'll add a few pictures to make it a bit more visual. Here is (part of) my input in the Weighted Overlay tool:

And this is my output. Notice how there are only classes 1 to 3, although my input rasters were classified into 1 to 4. Is this normal?


texts: ASCII encoding error when updating field content in SQL Server table

texts:
I know there are plenty of posts and blog explanations about unicode errors, but I still can't figure out to handle it in my particular case. So here is my problem: I am writing a Python script to  update records in an SQL Server table, with field content from a shapefile that has been edited in ArcPad.

I use an arcpy.da.UpdateCursor to update existing records.

I've put # -*- coding: cp1252 -*- on top of my script as I work with French characters.
Just doing so, I get the ASCII encoding error:

  'ascii' codec can't encode characters in position 0-1: ordinal not in
  range(128)

I've tried putting u before the field content:

...
with arcpy.da.UpdateCursor(DBtable, DBFields, where_clause) as DBCur:
    for DBrow in DBCur:
        ...
        DBrow[8] = u"{} - {}: {}".format(DBrow[8], date , AProw[8]) # AProw comes from a SearchCursor reading the shapefile.

I can then print a message with the text content, and it doesn't return an error message, but all characters are replaced (I get something like ???4???5????????>???5??? in my database field).

Any accentuated or punctuation character can be present in the fields (there's a free text comment field) so I don't want to check for every possible non-ASCII character and replace it.

I work with ArcGIS 10.2.2, the shapefile is edited in ArcPad 10.2 and the database is SQL Server 2008 R2.

What am I missing?

EDIT: This only occurs if the destination database is SQL Server. No problem with a file gdb. I have to add that the SQL Server table already contains non-ASCII characters.


texts: Warning displayed even though the file is detected in QgsMessageBar

texts:
I try to display a warning message when a file is not found in a folder but I don't understand why the warning message appears even if the layer is loaded.
I used if then else but it is not working.
from qgis.core import *
from qgis.utils import iface
import os

folder = "C:/Users/Desktop/working_folder"

os.chdir(folder)
wholelist = os.listdir(os.getcwd())

for file in wholelist:
    if "shapefile1.shp" in file:
            filefolder= folder +'\\'+file
            shapefile1 = QgsVectorLayer(filefolder,file[:-4],"ogr")
            QgsProject.instance().addMapLayer(shapefile1)
            print(filefolder)
    else:
        iface.messageBar().pushMessage("Error", "shapefile1 has not been found in :" + folder, level=Qgis.Warning)


texts: Decompose the polygon with polyline inside

texts:
Here is my question:

I have a polygon shapefile represent the whole country area.
A polyline shapefile(with a bunch of line boundaries) represent the province boundaries

So, how to decompose the polygon by the polyline inside and get the shapefile file with each province  independently?

I'm a new user of QGIS.

I have done the work of combine inner/outer boundaries using MMQGIS->merge layershows like this:

http://i4.tietuku.com/cfcd52e6886109fa.png

But the next step Change this polylines into polygons  I can't cope with

The post about cut the shapefile using digitizing tools didn't work for me
When I use vector -> Geometry tools -> polyline to polygon, the result shows like this:

http://i4.tietuku.com/7bcda0980879ff24.png

It seems that the inner boundaries had none effect.


texts: How to use the output file of a script as an input file for the rest of the script in arcpy?

texts:
If my question does not make sense or needs clarification please let me know. I am struggling to solve it so I was struggling to phrase it into a question.

I have this code here:

import arcpy,os
from arcpy import env

#Name .gdb and set location. This will create the .gdb
out_folder_path = ("C:/Users..../Documents/test")
out_name= "Test.gdb"
arcpy.CreateFileGDB_management(out_folder_path, out_name)

# This will load the .xml into the .gdb created above
target_gdb = "C:/Users/..../Documents/test/Test.gdb"
in_file = "C:/Users..../MASTER.xml"
import_type = "DATA"
config_keyword = "DEFAULTS"

# Execute ImportXMLWorkspaceDocument
arcpy.ImportXMLWorkspaceDocument_management(target_gdb, in_file, import_type)

I am struggling to convert this into a tool. The problem I am encountering is I want the user to just set the .env settings, then input the database name and then location. Then the tool will take over and create the file gdb and then load the .xml workspace into it. The .xml workspace has been hardcoded to a certain location. I am just struggling with the target_gdb for the .xml import. It does work as is with hardcoded locations and for one time use. But for those who do are not python saavy having to change code everytime is not feasible. Again I want this to become a tool.

I did find this formum: Automating script to set geodatabase?

but this just confused me even more.

Updated code:02/4/14

import arcpy,os
from arcpy import env

Name .gdb and set location. This will create the .gdb

out_folder_path = arcpy.GetParameterAsText(0)
out_name= arcpy.GetParameterAsText(1)
arcpy.CreateFileGDB_management(out_folder_path,out_name)

This will load the .xml into the .gdb created above

target_gdb = out_folder_path + os.path.basename(out_name)
in_file = "C:/Users/.../MASTER.xml"
import_type = "DATA"
config_keyword = "DEFAULTS"

Execute ImportXMLWorkspaceDocument

arcpy.ImportXMLWorkspaceDocument_management(target_gdb, in_file, import_type)


texts: QGIS: Saga Channel network: how to prepare SRTM data

texts:
I am trying to use the Saga toolbox (2.2.3) in QGIS (2.18) to create a channel network from SRTM data using this tutorial 1.

It works, but the result I get when using the DEM that comes with the tutorial are much better than what I get when using SRTM DEMs. The channels are mostly not connected to each other, no matter what I try. Is there a way to fix this, or is it because of the SRTM data being not good enough for this analysis?

So far, I am using the following workflow:

Catchment area (recursive): Step: 1, Method: Deterministic, Convergence: 1,100000

2: Channel network: Initiation theshold: Greater than 100000. I have also tried Less than 100000 and other values, but that does not help.

Do I need to prepare the SRTM data first before this analysis works?


texts: Categorizing GPS data based on cumulative time using ArcGIS Desktop?

texts:
I have GPS data which includes “Lon”, “Lat”, and cumulative Time from 0 (sec) as fields.
I want to add an integer variable new field, “Trip,” which starts from zero, and is incremented by one after the respondent stopped at some location more than 180 seconds.
Hence the “Trip” field may be like (0, 0, 0, 1, 1, 1, 1, 2, 2, …).

I would like to do this on ArcGIS if possible, but could not find effective solutions.


texts: Old GDAL/OGR version with QGIS from ubuntugis in 14.10?

texts:
Installing QGIS 2.6 from ubuntugis installs an old GDAL/OGR (ogrinfo --version
GDAL 1.10.1, released 2013/08/26) - it should be 1.11 for support of the Open FileGDB connector in QGIS 2.6
So I don't have access on features in a FileGeodatabase.

How do I update to GDAL 1.11 without conflicts. The ppa for trusty (14.04) seems to contain the newer GDAL.


texts: Google Earth Engine: When mean pixel value is blank, take nearest pixel value

texts:
I have a code for calculating the mean value of a variable, and reducing it to get values for 5 regions of interest. The full earth engine code for this is here.
However, I'm getting blank/null values for some of the pixels during some time periods. In these cases, is it possible to take the value of the nearest pixel?
My code is below. Is there an add on to this I could use, before reducing it by region?
// CURRENT //

//select current band from image collection
var currentCollection = ee.ImageCollection('HYCOM/sea_water_velocity')
.select('velocity_u_10','velocity_v_10')// east
.filterDate(startDate, endDate)
.filterBounds(EEZ); //EEZ

print('currentCollection', currentCollection);

// current speed needs scaling
function currentscale_function(image){
  var currentScale = image.expression(
 "velocity_u_10*0.001" &&
 "velocity_v_10*0.001",
 {"velocity_u_10" : image,
  "velocity_v_10" : image,
 });
   return currentScale.set('system:time_start', image.get('system:time_start'))
              .set('system:time_end', image.get('system:time_end'));
 }

var currentScale = currentCollection.map(currentscale_function);

    // calculate current speed by c=sqrt(a2+b2) using northwards and eastwards speeds
function current_function(image){
 var currentSpeed = image.expression(
   "sqrt(u*u + v*v)",{
   "u":image.select('velocity_u_10'),
   "v":image.select('velocity_v_10')
 });
 return currentSpeed.set('system:time_start', image.get('system:time_start'))
                .set('system:time_end', image.get('system:time_end'))
                .rename('current'); // rename band to current
}

var currentSpeed = currentScale.map(current_function);

var CURRENTmean = currentSpeed.mean();


texts: Set initial extent programmatically Android

texts:
I need to dynamically set the extent when the map loads.  I have graphics that are being added to the map layer and want the extent to allow all the graphics to be shown.

Currently I am getting a min x and y value and max x and y value based off the points that are being used to add the graphics.

Then I try to use set extent with an envelope that contains the min and max values of x and y.

So am I missing something to set the bounds so that the map will show the correct portion of the map?


texts: Is there a Time Slider feature in MapInfo?

texts: Is there an Time Slider Feature in MapInfo?  The only conclusion I have found in my research is that you have to buy the MapInfo Profiler arm?


texts: Setting a point using GDAL WFS-T tinyows

texts:
I’m trying to set a point xyz, which represents longitude altitude latitude, but when I create the feature, I get error:

  Error 1 Can’t find transcation response.

The Geometry field is set by this in the database:

ALTER TABLE poiserver ALTER COLUMN geom TYPE GEOMETRY(POINTZ,4326) USING ST_Force3D(geom);

OGRLayer *poLayer;
    poLayer = m_wfs_t->GetLayerByName("tows:poiserver");

    OGRFeature *poFeature;

    poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
    poFeature->SetField( "name", "ahmed" );

    OGRPoint pt;
    pt.setX( 120 );
    pt.setY( 240 );
    pt.setZ(10);

    poFeature->SetGeometry( &pt );

    if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
    {
        std::cout <<"error creating a feature"<<endl;

    }

    OGRFeature::DestroyFeature( poFeature );


texts: How to check for Selection in Geoprocessing Widget WebApp builder?

texts:
I published a gp tool and used the out-of-box geoprocessing widget.
When a user selects a layer to be processed, the gp process will call the FeatureService(FS) x number of times to get all the features in the layer. There is about 1 million features in the FS. This takes around 2 mins. Very annoying for a user. This is of course only happens when a user forgets to make a selection. Which is the main issue here.
Things I've tried.

In my GP script I stop if there is over 10,000 features, used the getCount() method. Still requests all the features before it submits the job.

On the client side I tried to tried to stop it my not sending the job if there is more than 10,000 features. Still to late in the step. It still calls like 20 - 30 xhr in the network window. In the code below I added the  if (inputValues.Layers.features.length < 10000)  statement.
executeGP: function(){
this._clearLastResult();
html.addClass(this.exeNode, 'jimu-state-disabled');
this._getInputParamValues().then(lang.hitch(this, function(inputValues){
this._showLoading();
  esriRequest({
    url: this.config.taskUrl,
    content: {
      f: 'json'
    },
    handleAs : "json",
    callbackParamName:'callback'
  }).then(lang.hitch(this, function() {
    if(this.config.isSynchronous){
      this.gp.execute(inputValues);
      this.infoTextNode.innerHTML = this.nls.executing;
    }else{
      if (inputValues.Layers.features.length < 10000) {
        this.gp.submitJob(inputValues);
        console.log("Good")
        this.tab.selectTab(this.nls.output);
      } else {
        alert("Please make a Selection")
        html.removeClass(this.exeNode, 'jimu-state-disabled');
      }
    }

  }), lang.hitch(this, function() {
    html.removeClass(this.exeNode, 'jimu-state-disabled');
  }));
}), lang.hitch(this, function() {
  html.removeClass(this.exeNode, 'jimu-state-disabled');
}))f

}

Ideas. I need to enable the run button only when there is a selection made? How do I got about doing this.  I know I didn't give all the code. But if someone could help me point me to the correct direction or have any other ideas how to solve this issue that would be helpful. Thanks!
I found this in the widget.js
    useSelectionAsInput: function(featureSet, layer) {
  array.forEach(this.inputNodes, function(node){
    if (node.param.dataType === 'GPFeatureRecordSetLayer' &&
        node.param.defaultValue &&
        utils.getTypeByGeometryType(node.param.defaultValue.geometryType) ===
        featureSet.geometryType) {
      node.inputEditor.setFeatureSet(featureSet, layer);
    }
  }, this);

  this.tab.selectTab(this.nls.input);
}

Not sure how this plays with the executeGP function. This happens in WAB 2.13. When I'm on WAB 2.16 it stops the call immediately not sure how, but I need this to work on 2.13
The image below shows that it's making all these requests.


texts: Counting lines in Python script with code block

texts:
I am writing a Python script using ArcPy and the Field Calculator function (code snippet below) and getting the following error message:

ERROR 000989: Python syntax error: Parsing error SyntaxError: invalid syntax (line 1)

Is line 1 the first line of the code block or the line calling the code block?
netYield = """netYield(ISlimit,multiplier,PASSA,compZone):
x = ISlimit / multiplier
if PASSA == 0:
    if compZone == "Office/Commercial":
        if x < 16000:
            return x
        else:
            return 16000
    elif compZone == "Industrial":
        if x < 111000:
            return x
        else:
            return 111000
    else:
        return x
else:
    return x"""

calcField(HDC,"NetYieldSQFT","netYield(!IS_limit!,!multiplier!,!PA_SSA!,!Comp_Zone!)", "PYTHON",netYield)


texts: How to see joins (one-to-many) in the QGIS 2.18.14?

texts:
I made a join with two tables - one to many.

territories and roads. Meaning, one territory have many roads.

territories -> 98 territories/features
roads -> 529 roads/features

So, in the territories, properties, made a join:

But in the attribute table I still see only 98 territories. Shouldn't see more than 98? Shouldn't see another rows for the related roads...?


texts:
I'm new to Leaflet and anything web related. I have got the following map working which asks for your location and changes the map view to that location. That is great but then the map continues to snap back to that location once you have panned round.
Is there a way I can prevent this?
<!DOCTYPE html>

    <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no" />
    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.3.1/leaflet.css" />

    <link rel="stylesheet" href="../dist/leaflet-gesture-handling.css" type="text/css">

© MapTiler, Contains OS data © Crown copyright and database right 2019


texts: ArcGIS 10.1 SA Raster calculator: Multi-raster conditional with nodata

texts:
I have a question regarding the Con statement when used in the SA Raster Calculator.

I am reclassifying a raster (categorical, integer) based on my own mappings but also on another raster (for some of the classes). So example I have

Con("pft.tif"==0,17,
Con("pft.tif"==1,1,
Con("pft.tif"==2,Con("aez.tif"==1,6,5),
Con("pft.tif"==3,3,
Con("pft.tif"==4,4,255)))))

I've written it like that for ease of reading.

My issue arises when there's no value for aez.tif (i.e. aez.tif value is NoData or the raster does not extend to the pft.tif extent).
From what I've seen the result is NoData. I thought it would result in "5". What I want is that if the Con with aez.tif results in this no value, that it does infact return some default that I've specified. I can't do that with e.g.

Con("pft.tif"==2,Con("aez.tif"==1,4,5),5)

since that terminates the Con while I still have other conditionals to try.


texts: Find Polygon in “Layer A” containing Polygon in “Layer B”

texts:
I haven't found yet someone with the same problem, so I'm posting a new question. The topic seems quite standard, so I'm sure there's a way to do it in QGIS.

I have two layers:

Region Layer Postcode Layer

both containing multi-polygon features.

I would like to cycle through all the Postcodes, find which is the Region that contains that Postcode, and eventually write the ID of the region inside the Postcode.

So far I've tried to write something in the Python Console, which I can share here:

from qgis.utils import iface
from PyQt4.QtCore import QVariant
_ID_FIELD = 'RS'
_NAME_FIELD = 'GEN'
_PLZ_FIELD = 'plz'
_NEW_ID_FIELD = 'PARENT_ID'
_NEW_NAME_FIELD = 'PARENT_NAME'
Land_layer = iface.mapCanvas().layer(1)
PLZ_layer = iface.mapCanvas().layer(0)
# Create 2 new fields in the layer that will hold the list of neighbors and sum
# of the chosen field.

#PLZ_layer.startEditing()
#PLZ_layer.dataProvider().addAttributes(
#        [QgsField(_NEW_ID_FIELD, QVariant.String),
#         QgsField(_NEW_NAME_FIELD, QVariant.Int)])
#PLZ_layer.updateFields()

## Create a dictionary of all features
PLZ_feature_dict = {f.id(): f for f in PLZ_layer.getFeatures()}
Land_feature_dict = {f.id(): f for f in Land_layer.getFeatures()}
# Build a spatial index
index = QgsSpatialIndex()
print 'PLZ Layer: %s' % PLZ_layer.name()
print 'Landkreis Layer: %s' % Land_layer.name()

for f in Land_feature_dict.values():
    index.insertFeature(f)
    print 'Inserting %s in Landkreis index' % f[_NAME_FIELD]

for f in PLZ_feature_dict.values():
    print 'Working on %s' % f[_PLZ_FIELD]
    geom = f.geometry()
    intersect_ids = index.intersects(geom.boundingBox())
    for id in intersect_ids:
        tempf = Land_feature_dict[id]
        print 'got parent Landkreis %s' % tempf[_NAME_FIELD]

Problem is, I'm getting the following error at line

intersect_ids = index.intersects(geom.boundingBox())

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "//VAIMUCO9/G114305$/Anwendungsdaten/Desktop/PLZ/ContainerFinder.py", line 34, in  <module>
    intersect_ids = index.intersects(geom.boundingBox())
AttributeError: 'NoneType' object has no attribute 'boundingBox'

Moreover, I'm not really sure that I should use these functions, and even this approach. Can anybody please suggest me some more clever way?


texts: How to query a WMS "singletile:true, ratio:1" in Openlayers 3/4?

texts:
I have a WMS that renders a logo into a corner of all requested images. I want the logo to be visible in the corner of the OpenLayers map.
I am upgrading from Openlayes 2.13 to Openlayers 4.6.4

In Openlayers 2 it worked by setting ratio:1, singleTile: true

  var layer = new OpenLayers.Layer.WMS(...,
  {
    singleTile: true,
    ratio: 1,
  });

In Openlayers 4 I don't get the "ratio:1" option to work.

How can I setup Openlayers 4 to request WMS-images with the exact same size of the map-view?

UPDATE: The current Single Image WMS example in OpenLayers 4 documentation shows a ratio: 1, but it doesn't seem to work for me.


texts:
I have to polygon shapefile maps in QGIS that cover the same area. One with ~90 polygons (call them regions) and the other with about ~800 polygons (call them zones). I want to add a data column in the zone shapefile that indicates what region the zone is located within. Most of the zones fall within a region, but zones fall on region borders. In this case, I'd like the variable to indicate the region that the majority of the zone falls in.

I tried searching around and I couldn't find any guides on how to do this. Any thoughts?


texts: Strategies for serving large datasets over WFS

texts:
We are implementing WFS using geoserver for our customers to consume vector data from a hosted spatial database to their desktop GIS.  They want to be able to just connect and browse the data in the same was that they would to for example an SDE connection.  Some of the datasets they want to consume are very large, and contain way more features than the feature limits we want to set on the server.  They may therefore not get all features downloaded, and if features are missing the user isn’t aware that they have hit the feature limit and don’t have all the features.

If we just increase feature limits, performance degrades, and the client still may not get all features, unless the limit is set larger than feature count in the data.

Basically we are experiencing performance and usability issues around guaranteeing that all features are downloaded for large datasets to desktop clients.

We are testing with Quantum, Arc9 and 10 with Data Interoperability connection.

So far we have been investigating how the caching options on the client can help.  We are coming to the conclusion that caching can help, but this still doesn't guarantee all features are always downloaded.

In quantum turning caching off and using scale thresholds helps, as features are requested again when the client zooms/pans if, but only if, the new envelope is outside the initial envelope.

In Arc 9 Data Interopeability, there is no ability to turn caching off, and no ability to restrict to an envelope, so the client gets the features on the first request for the whole dataset, limited by geoserver feature limit and never re-requests during the life of the cache(24hrs), so if features are missing, they will stay missing.

In Arc10 Data Interopeability, it's possible to use a live connection, which re-requests features on every pan/zoom, but as we increase feature limit on server, performance degrades.  Also, for certain tasks (view attribute table, select), arc10 requests all features one by one, which is pretty much un-usable.

In all cases, there are scenarios where the feature limit has been hit, features are missing, but the user isn’t aware because the client doesn’t flag this.

I would like some idea how people have approached serving large datasets over WFS to desktops in a way that guarantees the client gets all the features for their view extent.  Do we need to improve performance on the server and up the feature limit?  What are people doing to tackle this?

I would also like to try and get to a better understanding of what WFS is for.  Is WFS even appropriate for looking at large datasets in a desktop session as users are used to doing with shapefiles or a geo-database connections?

I realise this question is wide open, so any views much appreciated.


texts: How to set zoom level beyond 15 in Google Maps?

texts:
How do we set zoom level in Google Maps?

We can set zoom level to 15 but how can we set zoom level to more than 15?


texts: Import KML attributes into the correct schema

texts:
I have a kml file containing over 4000 polygons with attributes which must be retained (in 3 different fields). After using the standard KML to Layer tool, all the attributes are now in one PopupInfo field, along with the html code.

I have the Data Interoperability extension, so after looking at the questions here and here, I came across a post on the ArcGIS forums:

  The Data Interoperability extension will read KML (including folders and subfolders) and use any schema object defining the attributes. Thus what it writes out will have the correct field definitions. The same goes for writing KML with a schema object, given geodatabase inputs.

What I understand from that is that using the Quick Import tool will read the KML and export it to a gdb based on a schema which I defined. This seems to  be exactly what I want, except I cannot find the option to specify the schema.

I exposed the kml_name attribute as mentioned in the other question, but the field remained empty. I tried downloading the KML2shapefile plugin with MapWindow, but it is only for 32 bit and I am running 64 bit. I am unable to download the FME trial at work.


texts: including SRID in WKBWriter in Shapely doesn't work

texts:
I want to encode coordinates of LineString in WKB with the SRID included

 wkbwriter = shapely.geos.WKBWriter(shapely.geos.lgeos, big_endian=False, include_srid=True)

 geo = wkbwriter.write_hex(LineString(geometry_obj['coordinates']))

but it is the same exact as

shapely.geometry.shape(geometry_obj).wkb_hex

and does not include the SRID in the encoding, causing issues when trying to store it in postgis with column Geometry('LineString', 3857)

what am I doing wrong?


texts: Can you have a GetFeature response from GeoServer without coordinates in the response?


texts:
I have a remote sensing data in geotiff (input_data_rs.tif) format. I want to get corresponding ground truth from open source map (e.g. roads). I was able to download the mbtiles for the country. Is there a way to convert mbtiles to geotiff or geoJson (roads.geoJson) so I can use the following:

rio rasterize label.tif --like input_data_rs.tif < roads.geojson

I went through the https://github.com/developmentseed/skynet-data website but not quite sure how to use it for my data (input_data_rs.tif).


texts:
I'm using getGetFeatureInfoUrl to retrieve features from a WMS layer like that:

var url = wmsSource.getGetFeatureInfoUrl(
  evt.coordinate, viewResolution, 'EPSG:3857',
  {'INFO_FORMAT': 'text/xml'}
);

I have the problem that the features are quite small and so it's difficult to select a feature.

Is there any workaround to add something like hitTolerance in ol.interaction.Select? Or do I have to calculate manually an area and send a getGetFeatureInfoUrl request for every coordinate in this area?


texts: How to speed up QGIS project opening with large point datasets?

texts:
I recently loaded a big point set (in terms of memory usage) in the QGIS project i'm working on. This layer is necessary for the project but it is not used often. Most of the time is turned off.

The problem is that each time that I start the project, this layer is loaded, taking a long time for the project to start and making the project very slow.

Is there a way for not loading this layer into memory unless it is turned on?

I'm using version 2.12.3 but i guess it applies to older versions as well.


texts: Inserting attribute value in Print Composer into label based on spatial relationship

texts:
I have this print layout in QGIS print composer/atlas in which I only have one layer being shown: land parcels. I am interested in having a label that will show what is the attribute of another layer I have based on a spatial relationship like intersects.
This other layer is a layer of land uses. The idea is to have next to the map what are the allowed uses for that particular plot. The image below recreates what I need.

It seems there is no straight way to build this expression within the context of the QGIS Print Composer, but I think I am missing something. The closest I got an expression that feeds my label object is this expression:
if(intersects(@atlas_geometry, geometry(get_feature('land_use', "fid", "fid" > 0))),
  "use",
  'No restriction')

However this expression is always returning No restriction and this I am absolutely sure it cannot be. Therefore, this leads me to believe the intersects part is not being evaluated as I think it is. What could be wrong here?


texts:
I would like to add polygon geojson in the EPSG:2178 projection. The polygon should be displayed in central Europe but I see my geojson in Asia.
Where have I made a mistake? Where can I find information about ol.View "center"?

<!doctype html>

    <link rel='stylesheet' href='css/ol.css'>

    <title>OpenLayers 3 example</title>

My Map


texts: No valid training data were found [Google Earth Engine]

texts:
I would like to create a classification model through randomForest model. I used Sentinel 2 dataset. Everything is working if I'm using the original bands for the predicted variable. However, if I predict a custom variable (for example NDVI), I get this error:

  No valid training data were found.

I'm creating my model this way:

// Model variables
var output = 'B2';              // Working
var output = 'NDVI_GEE';        // Not working -> invalid data
// var output = 'NDVI_custom';     // Not working -> invalid data
var features = ['B4', 'B3'];

// Create random forest model + training
var model = ee.Classifier.randomForest({numberOfTrees:10})
                         .train(train,
                                output,
                                features);

// Apply classification
var predict = test.limit(50).classify(model)

Maybe is it a question of type ? I tried to cast them with no success.

See the full code on Google Earth Engine.


texts:
I am using GRASS GIS 7.6 on a CentOS server with 157 GB of RAM. The vector categories of the network layer that I need to run the Steiner tree algorithm on are as follows:

Layer/table: 1/network
type       count        min        max
point          0          0          0
line      364591          1     364591
boundary       0          0          0
centroid       0          0          0
area           0          0          0
face           0          0          0
kernel         0          0          0
all       364591          1     364591
Layer: 2
type       count        min        max
point       2174          1       2176
line           0          0          0
boundary       0          0          0
centroid       0          0          0
area           0          0          0
face           0          0          0
kernel         0          0          0
all         2174          1       2176

However even the 157 GB of RAM isn't sufficient and the process gets killed as follows:

GRASS 7.6.0 (srilanka):/home/asel/tmp/rds > v.net.steiner network arc_layer=1 node_layer=2 terminal_cats=0-2176 out=mysteinerSL
Number of terminals: 2172
Number of Steiner points set to 2170
Killed

I understand that this is due to the fact that the Steiner tree algorithm, at least at present, isn't very memory efficient. So I was wondering if I were to divide my network into segments and run v.net.steiner on each of them individually if I'd still get the desired output.


texts: Data source for determining postal codes given location?

texts: What is the best data source for determining international postal codes given a latitude & longitude?


texts: How to georectify a fictive earth to a projection-less version

texts:
I have a picture of a fictive earth-like world. Now I would like to draw the outlines and then georectify this to a projection-less version. I really need the map unrectified. Because I want to rectify it on my own.

Is that possible and even in QGIS?


texts: change width of outline for leaflet map depending on zoomlevel

texts:
I want to change the outline of a polygon depending on the zoomlevel. I thought it's easy to implement but there must be something wrong. My polygons dissapear when I try it with this code:

to get current zoomlevel:

map.on('zoomend', function() {
  var currentZoom = map.getZoom();
});

style function to style polygon:

function myStyle(feature) {
    return {
        fillColor: '#1c9099',
        color: 'white',
        if (currentZoom == 15) {
            weight: 2
        } else {
            weight: 3
        }
    };
}


texts: I want to do clustering of point data say buildings in an area to a fixed size for example each cluster should have 32x buildings. Can i do it in ArcGIS?


texts: SQL MAX Subquery in arcpy.MakeFeatureLayer_management

texts:
I am having an issue with subsequent subqueries in making a feature layer in arcpy.  My table has fields: OBJECTID, pointid (long), grid_code (double)

I can successfully extract MAX(grid_code), however it returns multiple values and as such I only need the row from selection where MIN(pointid) exists.

I have not been successful in using MakeFeatureLayer and then SelectByAttribute after that as arcpy effectively concatenates the two queries with AND, where I need to select using IN (I believe).

My code as currently stands below:

arcpy.MakeFeatureLayer_management(outpoint, "centroidpoint", "grid_code IN (SELECT MAX(grid_code) FROM point WHERE (SELECT MIN (pointid) FROM point))")

outpoint is a point feature class in a geodatabase that is the result of a raster to point conversion.

Error is:

  ExecuteError: ERROR 000358: Invalid expression grid_code IN (SELECT MAX(grid_code) FROM point WHERE (SELECT MAX (pointid) FROM point))
  Failed to execute (MakeFeatureLayer).

I believe error is in query, however I haven't found if I can do nested / subqueries of Max/Min funcitons.


texts: loading CanVec directory full of .shp files

texts:
I'm totally new to QGIS, and trying to view a Canadian Topo map as published on an iMac - I don't want to change it.

http://www.geogratis.ca/geogratis/en/option/select.do?id=5460AA9D-54CD-8349-C95E-1A4D03172FDF&datasets=031g05

the directory contains a lot of files - can all be loaded in one go so I get the whole map, not just pieces of it?
the CanVec attributes (colour, symbol etc.) aren't retrieved when loading a shape file with vector load, so that doesn't work to see it the way the map maker intended.

I suspect these problems are related, in that there's a file somewhere that lets QGIS tie everything together, but I can't figure out what it is or how to tell QGIS how to use it.


texts: Why is mapserver layer not showing up in openlayers?

texts:
I'm currently setting up a Mapserver/OpenLayers application. My mapscript works fine when i run on mapserver, but my mapserver layer not showing up in openlayers.

this is my mapscript "tangsel.map"

MAP
  NAME "TANGSEL"
  STATUS ON
  EXTENT 106.636718 -6.363274 106.779904 -6.228827
  SIZE 885 600
  SHAPEPATH "../shp"
  IMAGETYPE PNG
  IMAGECOLOR 252 253 255
  UNITS DD

  PROJECTION
    "init=epsg:4326"
  END

  LAYER
      PROJECTION
        "init=epsg:32748"
      END
      NAME "tangsel"
      DATA "tangsel"
      TYPE POLYGON
      STATUS ON
      CLASS
          NAME "tangsel"
          COLOR 0 0 120
          OUTLINECOLOR 0 0 0
      END
  END

END

file html "index.html"

            var map;

            function init() {
                map = new OpenLayers.Map({
                    div: "map",
                    projection: "EPSG:4326",
                    displayProjection: "EPSG:4326",
                    controls: [
                        new OpenLayers.Control.Navigation(),
                        new OpenLayers.Control.PanZoomBar(),
                        new OpenLayers.Control.ScaleLine(),
                        new OpenLayers.Control.MousePosition(),
                        new OpenLayers.Control.KeyboardDefaults()
                    ]

                });

                // create Google Mercator layer
                var gphy = new OpenLayers.Layer.Google(
                    "Google Physical",
                    {type: google.maps.MapTypeId.TERRAIN}
                );

                // create Mapserver layer
                var msv = new OpenLayers.Layer.MapServer( "SHP",
                    "http://localhost/cgi-bin/mapserv.exe",
                    {map: 'e:/ms4w/ppgt/map/tangsel.map'},
                    {layers: 'tangsel', srs: 'EPSG:4326', isBaseLayer: false, visibility: true},
                    {gutter: 15}
                );

                // Key Map
                var keymap = new OpenLayers.Control.OverviewMap({
                    maximized: true,
                    maximizeTitle: 'Show the overview map',
                    minimizeTitle: 'Hide the overview map'
                });
                map.addControl(keymap);

                map.addLayers([gphy, msv]);
                map.setCenter(new OpenLayers.LonLat(106.716677, -6.295165).transform('EPSG:4326', 'EPSG:3857'), 11);
                map.addControl(new OpenLayers.Control.LayerSwitcher());
            }

does someone know what the problem is?


texts: Mapbox Directions Plugin - Get Route GeoJson

texts:
Goal: find the GeoJSON points in a web map that are closest to a directions route using turf.js.

I am using the Mapbox Directions plugin. I am following this example.

How do I work with the route that is the output of the directions request?

Is there are way to get the route line in GeoJSON format to perform typical Mapbox functions on it?


texts:
I have the road networks from osm and I created topology.

I know that dijkstra algorithm works, because I calculated the shortest path from two point layers.
Now I would like to calculate catchment area for my points (from first layer).

I checked, how works pgr_drivingDistance. The result is far from expected.

I tried these codes:

SELECT pgr_drivingDistance( 'select id, source, target, leng as cost from osm_roads'
, 221, 1000)
     FROM osm_roads_vertices_pgr

where 221 is number of nodes. This code executed for a very long time and I had to stop it.

CREATE TABLE distances AS (SELECT a.node AS id, a.agg_cost AS distance, b.the_geom
FROM pgr_drivingDistance(
    'SELECT gid AS id, source, target, leng as cost FROM public.osm_roads',
    221,
    1000
) a, osm_roads_vertices_pgr b WHERE a.node = b.id);

but in this case result had only one point (start point) and this is not true.

I attach two screenshots:

The first, result from PgAdmin.

The second from visualization in QGIS.

What am I doing wrong?
Do you have an idea, how can I iprove my code?

I'm using PostgreSQL 10 with PostGIS 2.5.2 and Pgrouting 2.6.2.


texts: Reading geometry of a .gdb with large number of polygons

texts:
I've got a .gdb loaded into my ArcMap 10 that has about 9 Million polygons. I would like to get data from each polygon to include its geometry and certain fields from its row in the attribute table. I have done some reading here for geometry, and here for accessing the attribute table data, but I am unfamiliar with using the arcpy.da module along with cursors.

Could somebody provide me with an example of how I may loop through each polygon in the layer (it is the only one in my current map) and extract its linear ring of points and several field values?


texts: Creating opacity mask in Google Earth Engine?

texts:
I have an NDSI snow layer loaded in Google Earth Engine. Since snow is present for values <0.4, I want to make values >0.4 opaque. I looked at the manual's examples, but this code is not giving me what I want. Instead I get black (<0.40) and white (>0.40).

What am I missing?

var winter = NDSI.filter(ee.Filter.calendarRange(1999,2018,'year')).filter(ee.Filter.calendarRange(6,8,'month'));
var medianWinter = winter.median();
var medianWinterMask = medianWinter.gt(0.4);
Map.addLayer(medianWinterMask);


texts: SQL statement with GeoJSON driver in GDAL 2.2 not working

texts:
In GDAL 2.1.2 I was able to query from a GeoJSON data source using the following:

ogr2ogr -f geojson -dialect sqlite -sql "select * from OGRGeoJSON where NAME in ('Albania','Algeria')" countries.json world.json

Whereas in GDAL 2.2.4 this results in the following error:

ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(select * from OGRGeoJSON where NAME in ('Albania','Algeria')):
  no such table: OGRGeoJSON

What's changed?


texts: Why do I get a memory allocation error when selecting rasters by ID in PostGIS?

texts:
The Problem

I have about 256000 tiles in a table of GDEM v2 at 1 arcsecond resolution. I'm on an Ubuntu server running Postgres 9.5.5 and PostGIS 2.2. The table has a GIST index on the convex hull of each tile, and according to the docs, the primary key has a unique index created automatically so rid should be indexed as well.

I want to select a not-so-random subset of them for use in ST_Quantile.

The Query

I was successfully able to accomplish this using the following query:

SELECT (stats).*
FROM (
  SELECT ST_Quantile(
      ST_Union(rast),
      1,
      ARRAY[.1, .2, .3, .4, .5, .6, .7, .8, .9, 1]
    ) As stats
    FROM gdem_100

    -- this WHERE clause doesn't work
    -- WHERE rid BETWEEN 0 AND 500 OR
    --       rid BETWEEN 75000 AND 75500 OR
    --       rid BETWEEN 150500 AND 151000 OR
    --       rid BETWEEN 250000 AND 250500

    -- this WHERE clause works
    WHERE rid BETWEEN 0 AND 500 OR
          rid BETWEEN 1500 AND 2000 OR
          rid BETWEEN 3500 AND 4000 OR
          rid BETWEEN 4500 AND 5000 OR
          rid BETWEEN 5500 AND 6000

) As foo;

When I use the first (commented out) WHERE clause, I get a memory allocation error. The second one works fine.

One suggestion in the comments points to nodata values as a possible contributing factor ST_Union doesn't consider those by default. I will investigate this, but given the dataset, I think that quanitity of nodata would be really surprising.

Why does using the first WHERE clause cause a memory problem, when the second one, which actually selects more tiles, does not?


texts: How to run QGIS plugin from python script

texts:
I would like to access external QGIS plugins through a python script. I have been able to access the built in QGIS processing and vector toolbox, but have been unsuccessful with external plugins such as the topology checker plugin. I have tried this both using the built in QGIS python console, and an external IDE, but attempts have failed.

I am sure that there is a way to do this, has someone done this before?


texts: Receiving data from php script for geoJSON TileLayer?


texts: sf appears to truncate a projection string

texts:
Curious how folks deal with this type of thing where sf appears to trim the spatial projection:

library(sf)
library(sp)
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
st_crs(nc)$proj4string
[1] "+proj=longlat +datum=NAD27 +no_defs"

nc_sp <- as(nc, "Spatial")
proj4string(nc_sp)
[1] "+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat"

I am wanting to check if two projections are the same and I am running into this type of sf behaviour where the projections are the same but the strings aren't . Obviously I can deal with this via string manipulation but I am wondering if there is a more robust solution when projection strings are slightly different because of this type of truncation.


texts: Crash dump issue while uploading .csv coordinates as vector in QGIS

texts:
I am using QGIS 2.18 in Windows.

I want to show some coordinates in a map, for which I am uploading a .csv file, but it shows a crash dump every time with error "CRS was not defined". I have attached the image. Is there any solution please?


texts: Problems in changing projection using georeferencer in QGIS

texts:
I'm an amateur user and primarily use QGIS to change the projections of my fantasy maps, using the Georeferencer functon. Normally I enter non-equirectangular projections, plot the co-ordinates and export the image in an equirectangular projection. this is the extent of my experience using QGIS.

Anyway. I'm trying to do the opposite - I'm drawing an equirectangular map and have imported it into the QGIS Georeferencer, entered the co-ordinates and am trying to export in a non-equirectangular (Lambert equal area, if it matters), but it just keeps on exporting it as equirectangular. I'm not sure If I'm doing anything wrong. I've tried exporting in various different projections and after rendering, the exported image is always identical to the original one.

How can I make this work?


texts: If field 1 > field 2 update field 3 in loop

texts:
I am trying to update the table above.  What I want to do is if the values in the PMP fields are less than the Max fields update the check fields with PASS, if not FAIL.  I started doing this trying to create a list for the PMP fields and a list for the MAX fields and comparing them and updating a check list that I would use to populate the CHECK fields but had issues.

Could this be done using just an Update cursor or numpy somehow?

Here is what I had so far on my list attempt.

    pmpValues = []
    i = 0
    for pmpVal in pmpFields:
        with arcpy.da.SearchCursor(checkTable,pmpFields[i]) as cursor:
            for row in cursor:
                pmpValues.append([row])
        i += 1

    maxValues = []
    i = 0
    for maxVal in maxFields:
        with arcpy.da.SearchCursor(checkTable,maxFields[i]) as cursor:
            for row in cursor:
                maxValues.append([row])
        i += 1

    checkValues = []
    for p in pmpValues:
        for m in maxValues:
            if p < m:
                check = "FAIL"
            else:
                check = "PASS"
            checkValues.append(check)


texts:
I have been searching for an answer to this question for a while now, and for some reason I can't seem to find it.
How are the units defined in the EPSG:4326 system?
Are longitude and latitude written in degrees,minutes and seconds or are they written in decimal?
Some sources say that WGS84 and EPSG:4326 are exactly the same thing, other say they are different, but then don't provide the differences.
How should I be writing down longitude and latitude in the EPSG:4326 system?


texts: Converting CSV (polygon coordinates) to polygon feature class (shapefile)?

texts:
I have a CSV table as following :

ID;POLYGON
1;(-7.59443,33.55933):(-7.59416,33.55937):(-7.59414,33.55922)
2;(-7.58361,33.56655):(-7.58148,33.56655):(-7.58148,33.56486):(-7.58361,33.56486)

The coordinate system used is WGS84. Polygons have no holes and they are all single parts.

I need to convert it to a Polygon feature class/shapefile. I have access to QGIS or ArcGIS, but welcome answers in other software as well.

Any ideas on how to proceed?


texts: Hiding and re-arranging attribute table columns in QGIS?

texts:
This used to be possible by right-clicking column headings in the attribute table, then using the menu that appeared to produce a checklist of fields that could be unticked & re-arranged.

I've updated to 2.14.22 and the right-click functionality has vanished! Searched everywhere I can think of, can't find these settings or any clues in changelogs.

I'm sure I'm missing something here, can someone help me out?


texts: ArcGIS Pro not accepting negative volume value for Cut Fill

texts:
I am using the Cut Fill tool of ArcGIS Pro to calcuate volume change in DTMs over time. I am hoping to visualise these on a map where Cut (volume removed) is blue and Fill (volume added) is red. I was hoping to assign my own values for the primary symbology by classifying them in manual intervals. When I look at the attribute table, there is positive and negative values (as is to be expected). However, when I try and enter a negative value for the symbology classification, ArcGIS states that "The upper value specified is out of range. The default will be the minimum value." which is 1.0.
However, as you can see, there is both positive and negative values for VOLUME in the attribute table

But when I try and enter a negative value for the manual classification, ArcGIS prevents me from entereing a value below 1.0.


texts: 3D Raster's and Geosoft Montaj/Target, Leapfrog, GoCAD

texts:
I'm trying to figure out what the best format to serve 3D images inside our company.  Specifically georeferenced Tiff's that have been placed in 3d space with Geosoft Montaj.  We use ArcMAP extensively but don't want to go the way of the 3D Analyst extension.

We are beginning to store scanned Level Plans and Cross Sections and would like to store them in 3D, but what format?  The brass upstairs like using Leapfrog, but from what I can learn, it seems that Leapfrog doesn't like vertically geo-referenecd images.  It requires a manual entry of co-ordinates, elevations, and direction.  Which isn't exportable afterwards.

I've never used GoCAD before, but our geologists are evaluating it for geological modelling.  I've been able to look at it briefly and it's a bit intimidating.  Does someone have experience with it and importing 3d images?

So specfically, I'd like to know some people's thoughts who have used these programs for this purpose.  Also if there is a good 3d file format to get our data into so that it is as interchangeable as possible.  Thanks


texts: Reprojecting raster from one CRS to another?

texts:
I have a .tif raster file which is based on some projected coordinate system (for example: NAD 1983 State Plane California).
I would like to export this raster as a .xyz file, but also reprojecting it to a different projected coordinate system before export (for example to: UTM).

Do I have to reproject the initial .tif raster to WGS84 geographic coordinate system first, then reproject it to desired projected coordinate system (UTM), and then export it to .xyz file?

Or can I directly reproject it from one projected coordinate system to another one (UTM) and export it as .xyz file?

Essentially this is what I am would like to ask:

NAD 1983 State Plane California -> WGS84 -> UTM => export to .xyz

or

NAD 1983 State Plane California -> UTM => export to .xyz

I assume this problem can be applied to a general reprojection from one raster to another (without the need of exporting it to the .xyz file)?

I am using QGIS, but I assume the principle is the same in ArcGIS or other GIS software.


texts: Unable to change data source path for group layer file

texts:
I am having some difficulty changing the data source for a group layer. Below is the code that I am using. It processes a LYR file successfully that has a single data source but if LYR file has more than one data layer, it doesn't do it.

import arcpy, datetime, os, io, sys

try:
    #Read input parameters from GP dialog

    # Prod Setting...
    folderPath = arcpy.GetParameterAsText(0)
    output = arcpy.GetParameterAsText(1)
    sourceString = r"P:\Temp" #
    replacmentString = r"Z:\Temp" #

    #Create an output file
    outFile = open(output, "w")

    #Loop through each "LYR" file
    count = 0
    for path, dires, files in os.walk(folderPath):
        #for filename in os.listdir(folderPath):
        for filename in files:
            #fullpath = os.path.join(folderPath, filename)
            fullpath = os.path.join(path, filename)

        if os.path.isfile(fullpath):
            if filename.lower().endswith(".lyr"):
                groupLayerName = ""
                #Reference LYR
                layerFile = arcpy.mapping.Layer(fullpath)

                #Reference each layer in a layer file
                count = 1
                for lyr in arcpy.mapping.ListLayers(layerFile):
                    if lyr.isGroupLayer:
                        groupLayerName = lyr.name

                    if lyr.supports("DATASOURCE"):
                        dataSourceOrig = str(lyr.dataSource)
                        if dataSourceOrig.find(sourceString) >= 0:
                            dataSourceNew = dataSourceOrig.replace(sourceString, replacmentString)
                            if groupLayerName:
                                outFile.write(fullpath + "\t" + groupLayerName + "\t" +lyr.name + "\t" + dataSourceOrig + "\t" + dataSourceNew + "\t" + "GROUP LAYER SOURCE" "\n")
                                lyr.findAndReplaceWorkspacePath(sourceString, replacmentString, True)
                                lyr.save()
                            else:
                                outFile.write(fullpath + "\t" + "Non Group Layer" + "\t" + lyr.name + "\t" + dataSourceOrig + "\t" + dataSourceNew + "\t" + "SINGLE LAYER SOURCE" "\n")
                                lyr.findAndReplaceWorkspacePath(sourceString, replacmentString, True)
                                lyr.save()
                del layerFile

    outFile.close()

    #Open the resulting text file
    os.startfile(output)

    #Delete variables that reference data on disk
    del outFile

except Exception, e:
    import traceback
    map(arcpy.AddError, traceback.format_exc().split("\n"))
    arcpy.AddError(str(e))

(the code block above transcribed from this screenshot, indentation errors may remain)


texts: How to estimate the land surface temperature using ASTER image in ENVI 5.3?

texts:
I am doing research about land surface temperature using ASTER image.
I have read literature but i cannot do it in ENVI 5.3.
Could you provide step-by-step to calculate land surface temperature using ASTER in ENVI 5.3?


texts: What is the best way to convert WKT String to GML String

texts:
What is the best way to programmatically convert WKT String to GML String? Every reasonable library (at best available in some Maven repository) can be used.
Till now I came up with this solution using JTS. Its written in Java:
public String wktToGml(String wktGeometry) {
   WKTReader wktR = new WKTReader();
   Geometry geom = wktR.read(wktGeometry);

   // write JTS to string
   GMLWriter gmlW = new GMLWriter(true);
   String gml = gmlW.write(geom);
   return gml;
}

But I'm not sure if it's the correct and most effective way. There are probably some problems with LF (linefeed) symbols. Probably there should be also a solution using geotools library.
I'm new to this field, I'm trying to gain some overview.


texts: I am new using QGIS plugins. I used QtDesigner to build a plugin window with some buttons. But when I open this plugin in QGIS all the formatting is gone. Any idea what is happening?


texts: I have raster DEM and observer points (satellites) both are in the same coordinate system.The observer points are scaled down so that they are not very far away from the raster DEM (say at 500-800 m height). I want to do visbility analysis considering observer points as input and find out which are the locations in the raster surface am able i,e to see whether visible or not.I want to do it in arcgis or any other software


texts: Analysing local features on rasters using ArcGIS Pro?

texts:
I have several raster datasets representing the elevation on a beach (See Picture).
The beach topography include beach cusps (which represent local max on the beach).
I want to extract some statistical information from these features, including their height, their width, and their length in relation to their location on the beach.

I'm looking to find a way to extract this information using ArcGIS Pro but haven't been able to do so.

Does anyone has suggestions on how to do this in ArcGIS Pro?


texts: Converting non projected x,y coordinate into a specific projection in a 2D array

texts:
I created a 2D array and after applying Gaussian filter I find centroid of all the connected pixels with values equal to 1 and then
save it with the projection EPSG:3413.

Using skimage.measure.label, I determine all the connected pixels and label them. After labeling, I calculate their properties such as area and centroid using skimage.measure.regionprops().
I am able to extract the centroid coordinates of the connected pixels but these coordinates do not have any projection associated with it at this step.

Additional information:

How to add/associate projection to the 2D array, so that other attributes such as coordinates are in projection (EPSG:3413) and can be used for future analysis?

I know I can save the 2D array with projection (EPSG:3413) and then open and re-read the 2D array again and calculate coordinates which would have the projection. However I want to avoid this step of saving and then re-reading the 2D array as it would add more time to the computation.


texts: Using sql SELECT DISTINCT in Fusion Tables API

texts: How do I go about using sql SELECT DISTINCT in fusion tables API?


texts: Geocode address once in google maps, save in database, read geocoded result from there on

texts:
Is it possible to geocode an address, save the result in a database, then read that result from there on?

I'm using a Google map that is geocoding hundreds of addresses which obviously hit the geocode limit. Once geocoded, the address is very unlikely to change again so it makes sense to save the geocode result in a database.

I'm using Wordpress as a platform so if possible I'd like to save the geocoded result as a field in the database that relates to the particular company.

THE CLIENT HAS SPECIFICALLY ASKED FOR GOOGLE TO BE USED THROUGHOUT THE PROJECT THUS RULING OUT YAHOO OR ANYTHING ELSE.

UPDATE:
Everyone's saying about what I'm asking is going against their TOS but isn't their very own example doing the same thing - developers.google.com/maps/articles/phpsqlgeocode?


texts:
I have a linux machine which is connected by I2C to a F9P Ublock GNSS device. I need to get the GGA frames and some custom Ublocks frames from I2C, do some calculations and forward them to UDP.

Is there any simple way to use any standard service/library to achieve this parsing of NMEA/GGA frames or should I implemented my custom code for reading and parsing?

I have done some research and I am now trying to use GPSD.


texts: Performance problem with labels on query layer in ArcMap 10.7

texts:
I have a project in ArcMap 10.7.
In this project I defined a query layer that join a polygon feature class (SDE) with an Oracle view (non spatial).
For each polygon I have many (7) records in the joined view.
The query layer is time enabled.
All works fine if don't label polygons.
When labelling is activated redrawing becomes very slow (from 1 second without labels to 55 seconds with labels).
Labels are very simple, just one field.
Does anybody have any suggestions?


texts: Check if two different class points fall within the same cell of a grid?

texts:
I am very new to R and GIS and statistics. I'm trying to find out if two set of points of different class (ie. species) that I have are aggregated or seggregated. For that I need to randomize the proccess of the observation and resample to get a mean of how clustered are those observations.

First I created a grid over my study area and kept only the cells where observations of at least one type of points fell into and made sure the grids were squared polygons. That way I can check how many times the two type of points fall into the same cells compared to the total observed cells where points are found (Nº cells with coincidence /Nº of total cells with observations) and get a percentage. This was done in QGIS.

Then I need to randomize those observations and calculate again how many coincidences I have in relation to total cells with observations. If I repeat this process 1000 times I can get the mean of those 1000 result as the random aggregation pattern in that area to compare to the observed one. If my pattern was 60% and the observed 40% coincidences, then it would mean my points are more seggregated than randomly distributed in the area.

So, I want to do this in R. Firstly I uploaded the shapefile.

shape <- readOGR(dsn = ".", layer = "pruebamalla")

Then I generated random points within the cells of my shapefile for the first type of points. The number of points must be equal to my number of observations for that type of points, which are different for each type of points. (n is an example)

pointsinshape<-spsample(shape, n=120, type = "random")

And I can do the same with the other type of points sets

pointsinshape2<-spsample(shape, n=95, type = "random")

So now I have my two sets of points randomly distirbuted within the cells. Some cells would have x points of type 1, some x points of type 2, some x and y points of type 1 and 2, and some would be empty. I checked it with plot() and visually the results are reasonable

I know I can use over()from thesp package to check in which cells are my points of each type located. But now I don't know how to proceed to check which cells have the two type of points, besides checking manually in the list to see if the coordinates for the boundaries of each polygon given by over() match or not. And I need to know that to get how many coincidences I have, what is the total of different cells that have at least one type of point and to calculate the percentage of coincidences.

And if it can be done with a generic function or automatically so I could do 1000 resamples of that process and how to get the average valure of those results instead of doing it manually it would also be amazing as I need to repeat this process several times for different pair of points

I have already been searching everywhere and looking at similar threads or problems but found nothing, it would be amazing if someone could give me some guidance or orientation of how I could do this last step.

If you need more information or to upload the shapefile it isn't an issue.


texts: Getting data from Humminbird ChartSelect .bin file

texts:
I purchased some lake bathymetry data from humminbird chartselect. I can see those data using their applications such as FishSmart or also using the humminbird sonar module. But I would like to read those data as a CSV file with X, Y, Depth values. Is there any way to do so?

Here I'm attaching the screenshot of the data in their application.

Also, the file looks like this
Also, one of the bin file


texts:
I have a feature class that contains 131000 records. One of its fields (which was originally joined from a table) only has three possible values: null, 0, and 1. I have been trying to select the features that have a null value in this field, but each time I attempt to do so, it crashes. My syntax isn't incorrect ([FIELD] IS NULL). I'm not entirely sure the issue. Any ideas?

File Type: File Geodatabase Feature Class (field added from join from Excel table)
Version: ArcGIS 10.2.2 for Desktop


texts: How to import the "QGis" class in a script?

texts:
I want to check which type of geometry the current layer has. Therefore I use
the PyQGIS documentation on Access to Geometry and my code looks like:

if iface.mapCanvas().currentLayer().wkbType() == QGis.WKBPoint:
    print True

If I run this code I get following error:

Traceback (most recent call last):
File "<input>", line 1, in <module>
NameError: name 'Qgis' is not defined

I can't import a QGis class because I get an error that it does not exist.
Are there steps missing or how can I use this QGis class?


texts: Proper work flow to move manholes and pipes to surveyed location

texts:
I would like to know how to accomplish the following task using QGIS.

I have an existing point shapefile representing sanitary sewer manholes and an existing polyline shapefile representing sanitary sewer pipes.  All pipes are connected (snapped) to a manhole and all manholes are approximate location.  Recently a survey crew collected GPS points (4,000) for some of the manholes.  I need to move all manholes and pipes to their correct location.  How can this easily be done using QGIS?


texts: Deleting large temp files created in Image Analysis

texts:
I have a question about large temp files generated using Image Analysis in ArcMap.
I have been working with large digital elevation modesl and have been using Image Analysis to create hillshade and slopeshade maps for many (>2 GB) datasets. The files are supposedly stored as a temp file under appData and are equivalent in size. My C drive is filing up fast and I can't access the temporary folders. Is there any way I can get rid of these files so that I can free up disk space?
The path is: *AppData\Local\Temp\1\arcD660: but there are no files visible


texts: QGIS and remote PostgreSQL/PostGIS slow performance

texts:
I'm using a cloud service for hosting my PostgreSQL/PostGIS database.

When connecting with QGIS, its extremely slow.

The observed behavior is:

Connecting to PostGIS DB: works great.
Adding a layer to the map view: A bit slow but noting to worry
Open attribute table: really really slow, up to 45 seconds to open a small table of 100 point (QGIS become freeze and won't respond until the table open)

When Editing a layer in the PostGIS

Adding a new feature or editing the values of the attribute table: 45 seconds to get the form for insert the attributes and when saving the edit another 45sec. (sometimes can get worse if the geometry added its complex with alot of nodes its can take up to 3 min)

So I'm investigate a little and make some test, maybe the latency between my client QGIS (my PC) and the PostGis (cloud db) its the problem:

I'm using the Azure Database for PostgreSQL and the server is in US East region,
so I did the network test between client - host to see is the latency is the problem:

Using psql:

\timing
SELECT;
\watch 1

I get the following values:

83 - 84ms so it's not that bad, noting to get 45 sec in QGIS and the results are consistent with the ping from my location to the Azure server.

The DB performance looks ok since it's a really small db and queries with pgadmin o psql are fast.

Someone with experience using cloud db for using with QGIS maybe can guide me or give me an insight, maybe I'm missing a QGIS configuration or host server networking config.


texts: Resolving blanks when calculating attributes for almost 2 million polygons in QGIS?

texts:
I am starting with this Hazard Parcel Exposure dataset from the Metropolitan Transportation Commission which I downloaded as a shapefile.

http://opendata.mtc.ca.gov/datasets/parcel-hazard-exposure

I am using QGIS 3.4.3-Madeira on Mac OS High Sierra 10.13.6.

I tried to convert boolean to numeric values but the values did not calculate all the way through. Why are there blank rows when calculating attributes in QGIS3?

At first I thought maybe there is a selection on, but I did it all again with nothing selected and the same thing happened.

The 5 boolean attributes are shown below.

My goal is to convert each to a numeric value, sum into a new field, and represent each parcel with categorized symbology like so:

Where 0 means no hazards and 5 means there are 5 hazards (all five attributes are True).

First I added temporary fields for each hazard attribute and calculated so it is becomes numerical. This is what I calculated for the first boolean attribute.

"infltzone"='True'

The result of calculating that into the field "fltzone" is that True values become 1 and False values become 0. Then I sum all the values into a new field I called hazsum using this calculation:

 "fltzone"  +  "lndstdy"  +  "liqstdy"  +   "lndsus"  +  "liqsus"

The fields calculated and it initially looked like they calculated. But I tried to symbolize it and it looked like this, with many parcels without a numerical hazard exposure value:

When I identify a parcel with a hazsum, it looks great:

But when identifying a neighboring parcel with no hazsum...

I tried to export to Excel and do the calculations there but it kept crashing.

Can you suggest a method to accomplish my goal or to explain why my calculation did not complete?


texts: How to overlap and merge layers in QGIS?

texts:
I have two layers (1 and 2) and I want to overlap them and the intersection and the union didn't work.

I want to get something like this:

Is it really easy or is there a special method to follow ?

I am working with QGIS.


texts: How to stop the layer being removed from the map window after packing the table

texts: In MapInfo whenever I pack a table that's open in the map window, it closes it in the map window and I have to reopen it. How can I prevent this?


texts: How to merge multiple geotiffs by the maximum value using gdal_merge.py?

texts:
I'm try to use the gdal_merge.py to generate a a big composite image considering the maximum value per pixel (to avoid the no data value).

But I didn't find this option for gdal_merge.py.
My code is:

python gdal_merge -o mergedLayerOutput.tif -q -v --optfile filesToMerge.txt

The result:

Any suggestion about?


texts: Creating an Address Locator with a Database of Geo Information

texts:
I am trying to build an Address Locator for ArcGIS 10.1. I would like to use a Microsoft Server MySQL database where location information is stored using an alias, an address, and XY coordinates. What is the best way to create a geocoding service using this setup?

An option I have considered is writing a script to dump the SQL database into a format that plays well with ArcGIS on a regular basis but wanted to make sure there wasn't a better way to do this.


texts: Installing 64-bit geoprocessing of ArcGIS for Desktop?

texts:
I am trying to install 64-bit processing to ArcGIS 10.2, but the ESRI Patches and Service Packs website only has installations for Service pack 1 that goes with ArcGIS 10.1.

One similar question on this exchange shows a screenshot of an installation window that shows an option to include this, but my setup installation program (which I downloaded) does not show this installation window, but instead goes straight to selecting features to install as shown below.

I can't find any other documentation for running 64-bit geoprocessing on ArcGIS version 10.2.  Can anyone show me how to do this? or is it best to just try to install the service pack 1 version of it?


texts:
I am writing a standalone QGIS python application to run some unittests. I am getting an error when I try to use the QgsCoordinateReferenceSystem library.

ERROR 4: Unable to open EPSG support file gcs.csv. Try setting the GDAL_DATA
environment variable to point to the directory containing EPSG csv files.

I am running QGIS 3.0.3 on Windows 7.

I have made sure my GDAL_DATA path points to the gdal folder containing the gcs.csv. I have even placed a print statement above and below the call to QgsCoordinateReferenceSystem to print the environment variable which is correct as C:\Program Files\QGIS 3.0\share\gdal.

Any insight into this would be most helpful.

my full code is:

import os
from qgis.core import QgsApplication, QgsCoordinateReferenceSystem

app = QgsApplication([], True)
app.setPrefixPath(r"C:\Program Files\QGIS 3.0\apps\qgis", True)
app.initQgis()

app.setPluginPath(app.pluginPath().replace("\n", ""))
app.setPkgDataPath(app.pkgDataPath().replace("\n", ""))

print(os.environ["GDAL_DATA"])
c = QgsCoordinateReferenceSystem(4326)
print(os.environ["GDAL_DATA"])

This Python script is called from the following batch file:

@ECHO OFF

call "C:\Program Files\QGIS 3.0\bin\o4w_env.bat"
call "C:\Program Files\QGIS 3.0\bin\qt5_env.bat"
call "C:\Program Files\QGIS 3.0\bin\py3_env.bat"

set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%OSGEO4W_ROOT%\apps\qgis\python\plugins;%USERPROFILE%\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins%PYTHONPATH%

python test.py


texts: Creating building's shadow over buildings in QGIS

texts:
I want to create a map with the shadows projected from the building.

There are a few methods to create this, but all of them are not so accurate. The Translate Method has the problem that is not continous from the build. And the 2.5D Method is a quite good when I set this as a backward layer. But the problem is the shadow of a big tower is not over the small buildings, as well as the water tank box over the building not project the shadow.

I think this is a typical real 3D issue. But in QGIS is possible to solve this?

This image illustrates the question. The shadows should consider if the shape next is lower or higher.


texts: Removing a Table from several MXD-Files using ArcPy

texts:
I'm trying to remove a certain Table (xlsx-File) from some MXD-Files using python, but I can't finguring out how.

I tried the code sniplet for removing:

import arcpy, os
from arcpy import env

arcpy.env.workspace = 'V:/Projects/ZEL'

mxd =arcpy.mapping.MapDocument('V:/Projects/ZEL/test04.mxd')

df=arcpy.mapping.ListDataFrames(mxd, 'u')[0]

for lyr in arcpy.mapping.ListLayers(mxd, '', df):
    if lyr.name.lower()=='header_vorlage$':
        if lyr.dataSource == 'V:\Projects\ZEL\SHP\Daten_111004\ZEL_LB_Tabelle_V3_2010_2030.xlsx':
            arcpy.mapping.RemoveLayer(df, lyr)

mxd.saveACopy('pt.mxd')

print 'fertig'

del mxd

It works, but when I open the new created 'pt.mxd' there is still the symbol of the xlsx-Table with an '!' and you have the option to fix the datasource, because I have removed the path, but not the table...

I hope its not too confusing :)

My question is: Is there a way to remove a xlsx-Sheet from MXD-File without opening?


texts:
I am using QGIS version 3.2.2-Bonn on Windows 10 and this is my first attempt using Python on QGIS from the console. I am following this article:
Getting Started With Python Programming

When typing iface in the console a receive the following reply:

qgis._gui.QgisInterface object at 0x0000019EF2194D38

Then when I try to use the dir(layer) command I receive the following error message:

*dir(layer)*
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
NameError: name 'layer' is not defined

Do I need to configure something before starting using the console?


texts: Filling layer attributes field from Python list using PyQGIS

texts:
I have a QGIS3 layer which has an attribute column indicating the id of each features. I want to add two more fields and fill them with some values taken from a python list.

The list is something like that:

l = [['1', 'something', 'something else']]

where '1' is the id of the feature and 'something' and 'something else' are the values that must fill the column in correspondence with the right id.

    data_provider = layer.dataProvider()
    data_provider.addAttributes([QgsField('field1',  QVariant.String)])
    data_provider.addAttributes([QgsField('field2',  QVariant.String)])
    layer.updateFields()
    l = [['1', 'something', 'something else'],['2', 'something', 'something else']]

I tried in this way:

mem_layer.startEditing()
for feat in mem_layer.getFeatures():
    for i in l:
        if str(feat["id"])==i[0]:
            feat.setFields(feat.fields(),True)
            feat.setAttribute('field1', i[1])
            feat.setAttribute('field2', i[2])
mem_layer.commitChanges()

But nothing changed. How to solve this?


texts: Polygon rendering issue in map composer, QGIS 2.8

texts:
Since a picture is worth a thousand words, let's go straight to the point.
This is what I get in QGIS main window (and what I'd like to get in QGIS' map composer):

And this is what I actually get in the composer:

Note how the polygons representing the political claims (with colorized boundaries) carefully avoid the south pole in the map composer. What can I do? Googling the issue didn't help me much.

A few interesting precisions:

The same rendering issue sometimes arises, with this dataset, in QGIS main window. I can solve it by zooming in and then zooming out. This trick doesn't work in the map composer (I tried playing with the scale, with no luck).
When I downloaded this dataset, the "Simplify geometry" (Layer's properties > Rendering) was initially checked. Assuming this could be the cause of the issue, I unchecked it. The problem remains.

EDIT: the layer's CRS is EPSG:4326 (WGS84). It comes from Natural Earth.


texts:
When I do a timeline query using the Twitter API,

http://api.twitter.com/1/statuses/user_timeline.xml?screen_name=MyTwitterName

the geographic component of the reply looks like this:

<geo><georss:point>49.13786075 -122.85811069</georss:point></geo>

<coordinates><georss:point>49.13786075 -122.85811069</georss:point></coordinates>

Does anyone know the difference between the geo and coordinates tags? I did look in the Twitter docs but couldn't see anything. Should it always have both, or is one preferred to another?


texts: Resizing symbol changes its position on the map

texts:
I have character symbol for my point feature class from ESRI AMFM Water - Unicode 117 with offset Y:10,7 and size 36. At that moment it is placed perfectly as it is shown on a first printscreen. However when I change the size it gets shifted. I hope my printscreens represent perfectly what I mean. Do you how to change the way my points are resized?

I work with ESRI 9.3.1.


texts: Error: Algorithm not found with correct processing use

texts:
I try to run a simple dissolve algorithm in a standalone Python 3 script. My code is following :

from qgis.core import QgsVectorLayer

import processing
from processing.core.Processing import Processing
Processing.initialize()

layer = QgsVectorLayer(layerPath, 'layer 1', 'ogr')

parameters = {'INPUT': layer,
              'OUTPUT': 'memory:'}
dissolved = processing.run('gdal:dissolve', parameters)

But Python tells me : Error: Algorithm not found

But when I print the list of available algorithms :

print('\n'.join(map(str, [alg.id() for alg in QgsApplication.processingRegistry().algorithms()])))

I get the algorithm ID I try to use...

I run my scripts with a batch to include all libraries :

@echo off
set OSGEO4W_ROOT="C:\QGIS 3.0"
call %OSGEO4W_ROOT%\bin\o4w_env.bat
call %OSGEO4W_ROOT%\bin\qt5_env.bat
call %OSGEO4W_ROOT%\bin\py3_env.bat
@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%OSGEO4W_ROOT%\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%OSGEO4W_ROOT%\apps\qgis\python\plugins;%PYTHONPATH%
rem "%PYTHONHOME%\python" %*
%PYTHONHOME%\python .\main.py %*

Is there something I forgot to import or to do ?

EDIT :
after reinstalling QGIS I have another error message :

File "C:\Program Files\QGIS 3.2\apps\qgis\python\plugins\processing\tools\general.py", line 96, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:\Program Files\QGIS 3.2\apps\qgis\python\plugins\processing\core\Processing.py", line 125, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Error: Algorithm gdal:dissolve not found

But processing.algorithmHelp("gdal:dissolve") is working correctly...

Any idea ?


texts:
I have a shapefile with lines. In the attributes of the line shapefile I have known coordinates.

I want to get polygons out of these points which reach to the next line.
As a goal I want a grid out of polygons. Each polygon should contact each neighbour and should have the same ID with the known coordinate.


texts: Changing txt layer Data Source in QGIS only using Python console?

texts:
I want to change the data source of a layer (non geometric, just a table of fields) created from a .txt file. I know how to do it manually but I want to do it with the Python console. I have tried the setDataSource() command but I don't quite understand the arguments.

from qgis.core import QgsDataSourceUri

uri = QgsDataSourceUri()
uri.setDataSource("Path_to_new_source_txtfile","layer_name","memory")

This is the code I tried to use to change the source file, it didn't get me any error but it didn't change the source file either.


texts:
I have a map made of two layers, the first is a list of trails (path) and the second is a point based layer (different intersection between the trails).
I am looking for a way to select a segment between two point without selecting the all trail, here is a screenshot of what I'm having now, the yellow line is when I click, i get the whole trail instead of getting only the segment between two points (white dots), what I'm trying to do is having only the segment between the two point , where I click.
I am using mapbox gl js and turf.js and vector tiles.
also if you can help me find a better way of getting the second layer (intersection between the trails)

here is the code I made so far:
baseMap._map.on("click", "trails-data", (event) => {
  if (editMode) {
    var clickedFeature = baseMap._map.queryRenderedFeatures(event.point, {
      layers: [TrailLayerProperties.layerName],
    });
    var existingExtra = baseMap._map.queryRenderedFeatures(event.point,{
      layers: [SelectedLayerProperties.layerName],
    });

    if(existingExtra.length > 0){
      console.log(existingExtra);
      editComponent._selectedFeatures = editComponent._selectedFeatures.filter(
        (el) => el.id != existingExtra[0].id
      );

      baseMap._map
      .getSource(SelectedLayerProperties.sourceName)
      .setData(editComponent._clickedFeature);
      return;
    }

      let ark = [];
      let lines = [];
      editComponent._points.forEach((el) => {
        if (!turf.booleanDisjoint(clickedFeature[0], el.geometry)) {
          ark.push(el);
        }
      });
      console.log(ark);
      var ls = [];
      var feat;
      ark.forEach((it) => {
        var linesSplited = turf.lineSplit(clickedFeature[0], it);
        ls = [...ls,...linesSplited.features];
        console.log(ls);
      });
      let ptPt = turf.point(event.lngLat.toArray());
      if (ls.length > 0) {
        let min = ls[0];
        for (let i = 1; i < ls.length; i++) {
          let curMin = turf.pointToLineDistance(event.lngLat.toArray(),ls[i - 1]);
          console.log(curMin);
          let checMin = turf.pointToLineDistance(event.lngLat.toArray(),ls[i]);

          console.log(checMin);
          if (checMin < curMin) {
            min = ls[i];
          }
        }
        console.log(min);
        feat = min;
      }else{
        feat = clickedFeature[0];
      }
      feat.id = Math.floor(Math.random() * 100 * Date.now()); //random id generator

      editComponent._selectedFeatures.push(feat);
      editComponent._clickedFeature.features = editComponent._selectedFeatures;
      ark = [];
      lines = [];
      const length = turf.length(feat, { units: "miles" });

      baseMap._map
      .getSource(SelectedLayerProperties.sourceName)
      .setData(editComponent._clickedFeature);

  }
});


texts: QGIS 3.12 "Ballpark transform" shifted all the layers randomly on our project's due date!

texts:
I am using QGIS for data processing on a hydraulic model in another software that uses shapefiles for its data. The coordinate system is:
ESRI:102740 - NAD_1983_StatePlane_Texas_South_Central_FIPS_4204_Feet
and for some reason, QGIS always says "unknown CRS" when I open a file in this system, and I have to set CRS again myself. Today, I'm realizing that QGIS shifted over half of our shapefiles by a completely random transformation! This is what it has to say for itself:

An alternative, ballpark-only transform was used when transforming coordinates between Unknown CRS: BOUNDCRS[SOURCECRS[PROJCRS["NAD83 / Texas South Ce… and ESRI:102740 - NAD_1983_StatePlane_Texas_South_Central_FIPS_4204_Feet. The results may not match those obtained by using the preferred operation:
+proj=pipeline +step +proj=unitconvert +xy_in=us-ft +xy_out=m +step +inv +proj=lcc +lat_0=27.8333333333333 +lon_0=-99 +lat_1=28.3833333333333 +lat_2=30.2833333333333 +x_0=600000 +y_0=3999999.99999999 +ellps=GRS80 +step +proj=hgridshift +grids=ethpgn.gsb +step +proj=lcc +lat_0=27.8333333333333 +lon_0=-99 +lat_1=28.3833333333333 +lat_2=30.2833333333333 +x_0=600000 +y_0=3999999.99998984 +ellps=GRS80 +step +proj=unitconvert +xy_in=m +xy_out=us-ft
**Possibly an incorrect choice of operation was made for transformations between these reference systems. Check the Project Properties and ensure that the selected transform operations are applicable over the whole extent of the current project.**

How do I undo this "possibly incorrect" transformation?


texts: Using Formal Class RasterLayer as RSAGA input and otput

texts:
Is there a way to use the data from the R global environment object as input for an rsaga.geoprocessor alghoritm? And if yes, can be the output also saved as a same object?
I am writing chain processes and the constant importing/exporting of results makes the whole process far from elegant, as you can see on the code below:

rsaga.geoprocessor("grid_filter", 15, env = env,
                    list(INPUT="AKV_mode_5.tif",
                    OUTPUT="AKV_OUTPUT.sdat",
                    MODE="Moore",
                    THRESHOLD=4,
                    ALL="single class",
                    CLASS=2))

AKV_2_3 <- raster('AKV_OUTPUT.sdat')

AKV_2_3_recl <- reclassify(AKV_2_3, cbind(NA, 3))

writeRaster(AKV_2_3_recl, 'AKV_OUTPUT_2.tif', options="COMPRESS=LZW")

rsaga.geoprocessor("grid_filter", 15, env = env,
               list(INPUT="AKV_OUTPUT_2.tif",
                    OUTPUT="AKV_OUTPUT_3.sdat",
                    MODE="Moore",
                    THRESHOLD=4,
                    ALL="single class",
                    CLASS=3))


texts: Vector Tools Tab Does Not Appear on QGIS 3.10.8 and QGIS 3.14.1 on Windows 7 through Remote Access VPN (They do on Windows 10)

texts:
I am using QGIS on my work computer which runs Windows 7. I am using remote access (VPN) to access this work computer that runs Windows 7 since I am working from home. Today I tried to use the vector tools and I cannot access any of the options from the above ribbon (see attached photo). I have tried this in both QGIS 3.10.8 and QGIS 3.14.1 and in both instances, the vector tools don't show up when I press "Vector" in the above ribbon.
I have another computer at home that runs the latest QGIS on Windows 10 and the vector tool in the above ribbon shows the vector options without a problem.
I was wondering if someone has experienced this before? I assume this might be due to the operating system or VPN? What is strange is that last week the GIS Vector tools worked fine on Windows 7 that was accessed through remote access (VPN). Is there perhaps a plugin that is not downloaded?
I am more than happy to provide any more information needed to resolve this.


texts: Inconsistency with calculations of QGIS raster calculator and gdal_calc.py

texts:
I'm trying to compute a DEM using simple matrix algebra (DSM - DTM) but am running into a problem. The resulting .tif from the QGIS raster calculator is different to the output of gdal_calc.py.

The command I am using is:

gdal_calc.py -A nt6905_DSM_1M.tif -B nt6905_DTM_1M.tif --outfile=test.tif --calc="A-B" --format=GTiff --type=Float32 --debug

Following the answer in Why do I get different results with gdal_calc within the osgeo4w shell and QGIS' raster calculator?, I've changed the datatype parameter accordingly but am still returning different results.

Does anybody know how I can fix this? Is QGIS even using gdal_calc.py under the hood?

Sam


texts: How to get attributes of overlaying polygon layers by a spatial query?

texts:
I am new to QGIS and would like to merge selected data/attributes out of the attribute tables from overlaying vector objects on different vector layers without creating new additional  layers. Is there any standard QGIS operation or QGIS plugin to accomplish this task.

Basically to explain it better I have a vector layer with polygons representing buildings and its assigned attribute table. A second vector layer with polygons represents communal planning zones and its assigned attribute table. Both layers can not be joined by attributes as there have no corresponding keys. So I would like to perform a spatial query over both vector layers to find all buildings overlapping polygons of the zone layer and write back to the same building layer attribute table the corresponding zone code from the zone layer attribute table.

How can this be done? I am appreciating any hint.

Many thanks in advance


texts: Deciding which server best for geospatial data servicing?

texts:
I want to serve both vector and raster data sets of large volume (ex: for a country).

Right now I'm using GeoServer for vector and MapServer for raster.

Is there any literature or reference such that which server is better for servicng vector data and raster datasets?

I got some reference from http://anitagraser.com/2010/06/08/geoserver-vs-mapserver/ and Comparing different open source GIS servers?.


texts: VRP Multiple deliveries to one order, ArcGIS network analyst

texts:
I am working on a Vehicle Routing Problem (VRP), with network analyst in ArcGIS 10.5.1, for a number of depots and orders in the city Berlin (Germany). From the distribution centres (depots) trucks will supply supermarkets (orders).

The supermarkets have a specific number of tons they need to get delivered. This number is present in the attribute table. The number of tons differs between 11 and 106. The trucks have a capacity of 30 tons, which means that some supermarkets can be supply at once, but others need several truck loads.

My problem is that only supermarkets that ask for less than 30 tons are reached. All the supermarkets which require more than one delivery (more than 30 tons) are skipped by the trucks.
I already applied the route renewals to the trucks, so they are refilled after suppling a supermarket.

Does somebody know how I can let the truck also deliver to the supermarkets, which require more than 30 tons of goods?


texts: Exporting multiple shapefile layers to KML

texts:
In my code below, I am trying to export certain layers as KML files. It works, but it does not maintain layer properties such as color and thickness, nor does it export the attribute table associated with the layer. Can someone help explain the other parameters that QgsVectorFileWriter.writeAsVectorFormat takes so I can make this happen. I tried to read the documentation for it, but I couldn't find anything helpful.

from qgis.core import QgsVectorFileWriter, QgsVectorLayer
from PyQt4.QtCore import QSettings

layers = iface.legendInterface().layers()

for layer in layers:
    name=layer.name()
    if "Road" in name:
    output_layer = r"C:/Users/efrank/Desktop/Dump/" +  name + ".kml"
    QgsVectorFileWriter.writeAsVectorFormat(layer, output_layer, "utf-8", layer.crs(), "KML")

print "Layer to KML Conversion Complete"


texts: Projection conversion with Python/Pyproj: Canada Albers Equal Area to WGS

texts:
How do I get the geographical coordinates (lat and lon) from x,y values in Canada Albers Equal Area projection?

The range of x and y values as follows:

min(x) = -5171461.906673405
max(x) = 3372679.8809247203
min(y) = 209275.49942131247
max(y) = 5395930.584568907

I tried specifying EPSG code of Canada Equal Area as 102001 (see https://epsg.io/102001), but this throws an error in Pyproj.

RuntimeError: no options found in 'init' file

Gdalinfo for the tiff files of projected coordinate system and WGS are provided below. If you see the upper right corner value (0d38'48.11"W) in geographic coordinates of the projected coordinate system gdalinfo output, it looks like something is wrong -- totally off.

Coordinate System is:
PROJCS["Canada_Albers_Equal_Area_Conic",
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.2572221010042,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4269"]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["standard_parallel_1",50],
PARAMETER["standard_parallel_2",70],
PARAMETER["latitude_of_center",40],
PARAMETER["longitude_of_center",-96],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
Origin = (-5171461.906673405319452,5395930.584568906575441)
Pixel Size = (12034.002517743838325,-12034.002517743838325)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-5171461.907, 5395930.585) (165d 9'48.70"E, 46d33'32.99"N)
Lower Left  (-5171461.907,  197241.497) (145d24'33.39"W, 22d41' 7.36"N)
Upper Right ( 3384713.883, 5395930.585) (  0d38'48.11"W, 62d22'40.07"N)
Lower Right ( 3384713.883,  197241.497) ( 60d 4'49.15"W, 33d 7'48.19"N)
Center      ( -893374.012, 2796586.041) (114d45'33.46"W, 63d59'36.53"N)
Band 1 Block=711x2 Type=Float32, ColorInterp=Gray
NoData Value=-999

Gdalinfo of the WGS file is here

Driver: GTiff/GeoTIFF
Files: ABoVE.hydro.wy.OUT_PREC.1987.tif
Size is 476, 166
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-172.125000000000000,41.875000000000000)
Pixel Size = (0.249474789915966,0.248493975903614)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-172.1250000,  41.8750000) (172d 7'30.00"W, 41d52'30.00"N)
Lower Left  (-172.1250000,  83.1250000) (172d 7'30.00"W, 83d 7'30.00"N)
Upper Right ( -53.3750000,  41.8750000) ( 53d22'30.00"W, 41d52'30.00"N)
Lower Right ( -53.3750000,  83.1250000) ( 53d22'30.00"W, 83d 7'30.00"N)
Center      (-112.7500000,  62.5000000) (112d45' 0.00"W, 62d30' 0.00"N)
Band 1 Block=476x4 Type=Float32, ColorInterp=Gray
NoData Value=-9999


texts:
I am trying to join this data to this data (UTM32     Ebenen file =  VG250_GEM.shp) using columns Column2 and SDV_RS, respectively. I am using QGIS 3.4

I am following the instructions found here . When I try to perform the join, the fields in the csv file are NULL. Following the advice of DPSSpatial, I have used the refactor tool to change the data type (of the columns that I want to Join on) to integer.

Do the two files have to have the same number of rows?
Must all rows find a match in order for the Join to be successful?


texts:
I would like to apply a median filter to an image (.tif).
I looked for solutions in QGIS, and found a lot of possibilities to apply a mean filter, but none for a median filter.
So I decided to try something with the python console in QGIS, but I am quite new at this... The code I have for now is the Following one :

import numpy as np

from scipy import ndimage

import matplotlib.pyplot as plt

im_med = ndimage.median_filter(2018blackandwhite, 3)

2018blackandwhite being the image I want to apply the median filter on.

The error I am getting is "SyntxError : invalid syntax". I have tried to add the extension .tif after the name of the image, and I have the same error.

I am quite lost


texts: Let's have a list of SRS with proj4 definitions (e.g. spatial_ref_sys table). Is there a way to determine whether given SRS is geographic (lon,lat) or cartesian (x,y)?


texts: How do I get rid of that blue circle that is on the cursor when drawing a point in OpenLayers?

texts:
In OpenLayers, before placing a point, there is a small blue circle the moves with the cursor. I want to get rid of that and can't figure out how. If possible, it must be in the following section of the code:

var drawFeature = new OpenLayers.Control.DrawFeature(
  vectorLayer,
  OpenLayers.Handler.Point);
map.addControl(drawFeature);

I just can't seem to find anything on it!


texts: Listing Joined tables in MXD using Python?

texts:
I've a question about arcpy List* modules.

I have a quite large amount of data on ArcSDE Server.

Earlier I made a list of the feature classes, and tables in the mxd-s placed there using da.Walk and ListFeatureClasses/ListTableViews.

Now I'd like to list the joined tables too, but I've no further idea how to do it.

Actually is it possible to get these relations?


texts:
I have 45 maps produced with Atlas. I would like to save them with their names. However, there are character in these names that windows doesn't accept in file name. So, I'm tring to do a expression which replace this kind of characters to "_". I am doing it in output file name expression (panel in Genarate an Atlas).

So, here is my code. It isn't work what I want.

a = lower("field_name")
b = '-/"()' #characters to replace -, ", ( and ).
for char in a:
    a = a.replace(char, "_")
    'mapa_'||$feature||a