texts: Creating layer with WFS, but no features are added to the layer

texts:
I am creating a plugin that is gathering data from a WFS client. The WFS client was recently updated, and I had to update the URL and filter encoding that I was using. I updated the code, and it is now working fine on one computer, but when I later tested on another computer, it did not. It creates a layer as normal, but no features are put into it. I think my URL is not encoded correctly, but I am not sure what is incorrect in it.

I create the layer with QgsVectorLayer("url", "name, "ogr") and then QgsMapLayerRegistry.instance().addMapLayer(new_layer). I am trying to encode the URL by changing some of the values:

< : %3C
> : %3E
space : %20
" : %22

An example of a URL and filter encoding that works on one computer, but not others:

url: http://wfs.geonorge.no/skwms1/wfs.tilgjengelighet_tettsted?service=WFS&request=GetFeature&version=2.0.0&srsName=urn:ogc:def:crs:EPSG::3857&typeNames=app:TettstedInngangBygg&

filther encoding: FILTER=<fes:Filter xmlns:app="http://skjema.geonorge.no/SOSI/produktspesifikasjon/TilgjengelighetTettsted/1.2"><fes:PropertyIsEqualTo><fes:ValueReference>app:kommune</fes:ValueReference><fes:Literal>0301</fes:Literal></fes:PropertyIsEqualTo></fes:Filter>

After I have changed some of the characters to encode it, it looks like this:

http://wfs.geonorge.no/skwms1/wfs.tilgjengelighet_tettsted?service=WFS&request=GetFeature&version=2.0.0&srsName=urn:ogc:def:crs:EPSG::3857&typeNames=app:TettstedInngangBygg&FILTER=%3Cfes:Filter%20xmlns:app=%22http://skjema.geonorge.no/SOSI/produktspesifikasjon/TilgjengelighetTettsted/1.2%22%3E%3Cfes:PropertyIsEqualTo%3E%3Cfes:ValueReference%3Eapp:kommune%3C/fes:ValueReference%3E%3Cfes:Literal%3E0301%3C/fes:Literal%3E%3C/fes:PropertyIsEqualTo%3E%3C/fes:Filter%3E

If i run them in a web browser, I get results, but it also says in the top: numberMatched="unknown" numberReturned="0"

I am hoping someone can see what I am doing wrong or are able to point me in the right direction. I have lost a lot of time because it appears to be working on my main computer, and I didn't have a lot of time for testing.


texts: Converting GMT grid to ArcGIS readable format?

texts:
I am trying to convert a GMT grid data into a format that is readable in ArcGIS. I have already attempted to load the data into GRASS and use GDAL translate to convert into a tif with unwanted results. Mainly the extent is wrong with a cell size of 1,1 while in reality it should cover the entire world.

Update

After installing GMT and utilizing the source code

grd2xyz -Ef input.grd

the resulting ascii grid was converted to a tif format that was read by ArcGIS (ascii to raster). However again the final product does not have a correct extent showing a cellsize of 1,1. Could someone who perhaps is more familiar with GMT check the data above to check if it is viable to use within an ArcMap format?


texts:
I found OpenJump a real nice tool with its ability to render results of ad hoc SQL queries like PostGIS ones (among other things).

However, I also found that not everything that runs in (say) pgAdmin can be executed in OpenJump (File > Run Datastore Query). For example:

If I run a query that does not return result,

CREATE TABLE t AS
SELECT * FROM tbl;

I get an error (as of Version 1.10 release rev.5336) about No results were returned by the query. Invalid query .... (Yet if I run the above again, it says Error: relation t already exists. Seems that the query is executed anyways.) Similarly problem happens with CREATE FUNCTION etc. The error happens even if there is a second query that does return a relation.

If I run the same query twice,

SELECT * FROM tbl;
SELECT * FROM tbl;

I get an Error Multiple ResultSets were returned by the query. Invalid query....

From trial and error, it seems that the only things found to work so far is a single SELECT statement or equivalent (CTE etc.).

I am wondering if there is any documentation on this. In particular,

What kind of queries are valid?

How many queries can be performed in one run of data store query?

(Wouldn't it be nice if OpenJump just returns the result of the last query?)


texts: How to use Arcpy UpdateCursor to update a field?

texts:
This should be simple but I've spent days trying everything to get it to work. I need to update a field in a layer based on two other fields. I can't even get it to update based on one field. Could someone please tell me whats wrong. So in this script I want to be classifying the outfield as "poor" if the infield is "Sandy". I will also have another field of numbers (and infield2 = "10").

rows = arcpy.UpdateCursor(infc) #this is my feature layer
for row in rows:
    if row.getValue(Infield1) == "Sandy":
        row.setValue(Outfield, "poor")
    else:
        row.setValue(Outfield, "none")
    rows.updateRow(row)
del rows

Thanks for all the answers. The script only works if written like this: I did not want to name the field names though, but at least it works.

The only way I could do a range 5 - 10 was to put > 5 and <10 (see line 23)

for row in rows:
    if row.getValue("Texture") == "Sandy" and row.getValue("Humus") < 10:
        row.setValue("Richness", "poor")
    elif row.getValue("Texture") == "Sandy" and row.getValue("Humus") > 10:
        row.setValue("Richness", "medium")
    elif row.getValue("Texture") == "Sandy-Clay" and row.getValue("Humus") < 5:
        row.setValue("Richness", "poor")
    elif row.getValue("Texture") == "Sandy-Clay" and row.getValue("Humus")  >= 5 and row.getValue("Humus")  <= 10:
        row.setValue("Richness", "medium")
    elif row.getValue("Texture") == "Sandy-Clay" and row.getValue("Humus") > 10:
        row.setValue("Richness", "rich")
    elif row.getValue("Texture") == "Clay" and row.getValue("Humus") < 5:
        row.setValue("Richness", "poor")
    elif row.getValue("Texture") == "Clay" and row.getValue("Humus") > 5:
        row.setValue("Richness", "rich")
    else:
        row.setValue("Richness", "unclassified")
    rows.updateRow(row)
del rows, row


texts:
I am using ArcGIS 10.3.1. I have feature classes in a geodatabase. I am mapping noxious weed patches, and would like to auto-populate the ownership info, mailing address, etc.  The "parcels" feature class contains this info. Parcels may contain more than one patch. Each parcel has a unique "PIN" and I would like to auto-fill by entering the PIN into the patches attribute table.

Is this possible?  What is the correct tool for the job?


texts:
Here's my code:

        function main() {
            // create leaflet map and define some properties
            var map=    L.map(
                "map",
                {
                    zoomControl: true,
                    center: [40,-95],
                    zoom: 4
                }
            );

            // add a base layer to map
            var Hydda_Full = L.tileLayer('http://{s}.tile.openstreetmap.se/hydda/full/{z}/{x}/{y}.png', {
                attribution: 'Tiles courtesy of OpenStreetMap Sweden &mdash; Map data &copy; OpenStreetMap'
            });
            Hydda_Full.addTo(map);

            // viz.json of map
            var layerURL=           "https://CARTO_USERNAME.carto.com/api/v2/viz/CARTO_ID/viz.json";

            // add CARTO layer with on layer
            cartodb.createLayer(map, layerURL, {infowindow: false, tooltip: false})
                .addTo(map)
                .done(function(layer){
                    layer.setInteraction(true);

                    layer.on("featureClick",function(e, latlng, pos, data){
                            console.log(data);
                        });
                    });
        }

        window.onload = main;

console.log(data) prints an object with only the cartodb_id key and its value. But the dataset on which my Carto map is based has way more columns than cartodb_id.

Why do I only get the cartodb_id column, not the others?


texts: OpenLayers 6: add padding on fit function

texts:
I'm trying to apply the padding on fit function but is seems that there aren't effect. This is my code:

var sourceGeoBlog = new ol.source.Vector({
  url: '{% url 'api_geopost' %}',
  format: new ol.format.GeoJSON(),
});

var geoBlog = new ol.layer.Vector({
  source: sourceGeoBlog,
  style: styleFunction,
});

sourceGeoBlog.once('change', function(evt){
  if (sourceGeoBlog.getState() === 'ready') {
    if (geoBlog.getSource().getFeatures().length > 0) {
      var padding = [500, 0, 500, 0]
      map.getView().fit(
        sourceGeoBlog.getExtent(),
        map.getSize(),
        {
          padding: padding,
        }
      );
    }
  }
});

var view = new ol.View({
  center: ol.proj.fromLonLat([4.0863647111, 48.6698370737]),
  zoom: 9,
  minZoom: 7.5,
});

var map = new ol.Map({
  controls: ol.control.defaults().extend([
    new ol.control.ScaleLine({className: 'ol-scale-line', target: document.getElementById('scale-line')}),
    new ol.control.LayerSwitcher({tipLabel: 'layer switcher'}),
  ]),
  target: 'map',
  layers: [
    basemapGroup,
    geoBlog,
  ],
  view: view,
  overlays: [overlay],
});

I can zoom on extents of geoBlog.
What I'm doing wrong?


texts: QGIS 2.16.2 installation on Windows: 'Extract: error writing to file_flapack.pyd'

texts:
I encountered an error while installing QGIS 2.16.2 Nodebo (see below)
and I don't know how to address this problem.


texts: Calculating minimum bounding circle in QGIS or Python?

texts:
I'm attempting to calculate the radii required to cover a set of polygons in a shapefile; for example, given any number of neighborhood polygon boundaries, I'd like to calculate the smallest circle radius required to cover each.

Is there a way to calculate a minimum bounding circle in QGIS, given a polygon? Otherwise, is there a way to do so in python?

Edit: I'm using QGIS 2.18 on MacOS.


texts: How to Hide InfoWindow/Popup when "no information available" or when buffering?

texts:
I have an app that does some buffering, only problem is when I click on the map to put in the point/line/polygon vertices to do the buffer, I get popups at every mouseclick. It doesn't interrupt the buffer process but its super annoying.

At first thought I wanted to try and disable the popups so that they didn't show when "No information available" was the result. But I figure there will be more times than not when a popup will return some results and thus, need to disable/hide the popup so that it doesn't appear at all while the buffer process is running.

I've tried putting map.infoWindow.hide(); within the buffer without success and also saw this post recommending an onClick handler and have been unsuccessful as well. Any tips appreciated.


texts:
I want to zoom to ol3 map base on scale(for example 1:50000). This means every 1cm in ol map must show 50000cm in real world. for converting scale to resolution I used as follows:

var getResolutionFromScale = function(scale, units) {
    var dpi = getDpi();
    var mpu = ol.proj.METERS_PER_UNIT[units];
    var inchesPerMeter = 39.37;
    return parseFloat(scale) / (mpu * inchesPerMeter * dpi);
}

function getDPI()
{
    var div = document.createElement("div");
    div.style.width="1in";
    var body = document.getElementsByTagName("body")[0];
    body.appendChild(div);
    var dpi = document.defaultView.getComputedStyle(div, null).getPropertyValue('width');
    body.removeChild(div);
    return parseFloat(dpi);
}

Now for testing it, I use ScaleLine and when scaleLine controle show 1000m, it's length must be 2cm, but it is almost 2.5cm.

Where is th problem? How can I zoom base on scale?

Online Demo


texts:
I just wanted to gauge the work flow for adding charts to a map using the SDK that comes with OpenGeo suite. I know it seems a vague question I'm just trying to get my head around it :/

So far the documentation is pretty straight forward with adding plugins using the gxp.plugins such as  * @require plugins/FeatureEditor.js

However Im not sure about how I would go about adding geoEXT additions such as Link from OpenGeo Suite to GeoEXT charts. What would the basic starting point be, and any pointers would be great. I will eventually be looking to integrate into a popup on my map a graph, that would take values over a time frame and plot them.


texts: problem in PostGIS Name of the system used for Points and how to get distance between to point?

texts:
I have started to work on a project in PostGIS.

I have a table with a geometry points in which points saved like this :   "010100000089F511C671C01E4162327E2A83545341"

and I have translated this point in postgis using :
 ST_AsText(geometry g1); and the result was : "POINT(494852.884569681 5059777.30183276)"

I dont know the name of this systems.....I want to get the distance between to point ??

Would you please help me what is the name of this systems and how can I get the distance between points ?


texts: Obtaining correct area values *With* correct projections?

texts:
I'm trying to get the percentage of each U.S. county overlaying groundwater (data from the census and USGS respectively). They're both polygons.

I have 1) correctly projected the data and the layer (NAD 83, meters) 2) clipped the counties using a polygon corresponding to groundwater coverage 3) use calculate geography to obtain the areas, so that I can then figure out proportions/percentages. I create a new field for that (double).

I've been toying around with this and other methods, and I cannot get credible area values. They're tiny; almost all 0.01 - 0.1. One of my many hunches was that it's in decimal degrees, but I've projected the data, and it's still showing up incorrectly. The most frustrating part is that the county/state files won't calculate correctly either on their own, which is bizarre (Florida = 15 m2).

Any ideas?


texts: Get code to run on every column in a table -Error "Could not find a public static method named ..." -Mapbasic 17

texts:
I need to run a table update command on every column in a table. The code works fine when I manually define the column to work in but not when it is part of a for loop that gets all the column names I get

If I take the update table command out of the loop and into the code after the for loop, it works fine.

What's different with defining a declare for use in a loop?

'64 bit regex based on https://community.pitneybowes.com/communities/communityhome/digestviewer/viewthread?MessageKey=a92a0fcd-5348-4eae-aa93-a47d5a1e6bbe&CommunityKey=3c2aca7a-c3ae-4602-a142-9ee956769d55&tab=digestviewer#bma92a0fcd-5348-4eae-aa93-a47d5a1e6bbe'
Include "MapBasic.def"

'64 bit regex'
Declare Method RegExReplace Class "WindowHelper.Controller" Lib "WindowHelper.dll" (ByVal sInput As String, ByVal sPattern As String, ByVal sReplace As String) As String

Declare Method RegExIsMatch Class "WindowHelper.Controller" Lib "WindowHelper.dll" (ByVal sInput As String, ByVal sPattern As String) As Logical

Declare Method RegExFirstMatch Class "WindowHelper.Controller" Lib "WindowHelper.dll" (ByVal sInput As String, ByVal sPattern As String) As String

Declare Sub Main

Sub Main
    Dim inFile, regexString as String
    Dim inColumn as Alias

    inFile= "Z:\My Drive\Mangoesmapping\Spatial Projects\2018\DSC\007_DSC_Asset_Identification_and_Geocoding\Working\Phase3_DataEnrichment\Working\DataUpdates\Scratch\Water_Hydrants.TAB"

    'Perform Regex and text Replace'
    Open Table inFile as tempTable
    'Alter Table tempTable (add ScratchField Char(254) ) Interactive
    Browse * From tempTable
    regexString ="^\s\|\s(.*)"

    ' Perform on all fields -ref: https://gis.stackexchange.com/questions/75494/how-to-loop-through-all-columns-in-a-table-layer-getting-all-data-for-a-row-feat
    Dim sCol As String
    Dim aCol As Alias
    Dim nCol As Integer
    For nCol = 1 To TableInfo(tempTable, TAB_INFO_NCOLS)
        sCol = ColumnInfo(tempTable, "COL" & nCol, COL_INFO_NAME)
        aCol = tempTable & "." & sCol
    'Update every column as it's found --- GIVES ERROR
        'Update tempTable Set ScratchField = RegExReplace(aCol,regexString,"")
        Next
    'Update a known column
    inColumn = tempTable & ".As_con_3_hyperlink"
    Update tempTable Set ScratchField = RegExReplace(inColumn,regexString,"")
End Sub

UPDATED CODE based on T_Bacon's input.

Had to change scratchField to aCol as it is now checking every column.

Sub Main
    Dim inFile, inField, regexString, strValue as String

    inFile = "Z:\My Drive\Mangoesmapping\Spatial Projects\2018\DSC\007_DSC_Asset_Identification_and_Geocoding\Working\Phase3_DataEnrichment\Working\DataUpdates\Scratch\Water_Hydrants.TAB"
    inField = ".As_con_3_hyperlink"
    'Perform Regex and text Replace'
    Open Table inFile as tempTable
    'Alter Table tempTable (add ScratchField Char(254) ) Interactive
    Browse * From tempTable
    'Pattern for starting | without anything before it
    regexString ="(.*)[\s]+\|\s(.*)"

    ' Perform on all fields -ref: https://gis.stackexchange.com/questions/75494/how-to-loop-through-all-columns-in-a-table-layer-getting-all-data-for-a-row-feat
    Dim sCol As String
    Dim aCol As Alias
    Dim nCol As Integer
    For nCol = 1 To TableInfo(tempTable, TAB_INFO_NCOLS)
        sCol = ColumnInfo(tempTable, "COL" & nCol, COL_INFO_NAME)
        aCol = tempTable & "." & sCol
    'Update every column as it's found --- GIVES ERROR ref: https://gis.stackexchange.com/questions/316048/get-code-to-run-on-every-column-in-a-table-error-could-not-find-a-public-stati
        strValue = aCol
        Update tempTable Set aCol = RegExReplace(strValue, regexString, "")
        Next

End Sub


texts: Separating a certain section of lines in QGIS

texts:
I'm trying to separate a specific area of lines in a topo shapefile in Qgis so I can import it into Civil 3D without it taking several centuries to load in the full file, however when I use the "split features" tool it takes forever just to cut one line, and trying to cut multiple lines just generates an error.

Also, since it's hard to explain what exactly I'm trying to do in a few words: just imagine having a square polygon and cutting out a smaller square from the center and deleting the rest (but in this sentence replace "square" with "topo line data").


texts: How to fill spatialite geometry column with values?

texts:
This is a follow up to this question: Export vector layer to sqlite database using pyqgs

I am trying to import a QGIS polygon layer into a spatialite db.
I was able to add a geometry column to my table using

cur.execute("SELECT InitSpatialMetadata()")
cur.execute("SELECT AddGeometryColumn('Polygonx', 'geometry', 4326, 'POLYGON', 'XY')")

Subsequently, I am using the following line to insert my values, inlcuding the geometry

cur.executemany("INSERT INTO Polygonx VALUES (?,?,?,?)", feature_list)

Whereas

feature_list = []
for feature in vlyr.getFeatures():
    attributes = feature.attributes()
    new_feature = (feature.id(), attributes[0], attributes[1], feature.geometry().exportToWkt())
    feature_list.append(new_feature)

When using the db manager to check the table, however, the value of the geometry column is always NULL. Any idea what wrong with my code?


texts: Get pixel values from a raster using a multipartpolygon

texts:
I'm trying to get the mean value of pixels from an NDVI raster, that intersect the geometry of a FeatureCollection, that is a multipart polygon which represents my sample area, but I can't find the way.

Here's the code:

 // Load geometry
var limites_BSHC = ee.FeatureCollection ('users/charlieswall/tesis/BSHC_LIMITES');
var poligonos = ee.FeatureCollection('users/charlieswall/tesis/POLIGONO_BSCH_DATOS_COMPLETO');
var p_inicio = ee.FeatureCollection('users/charlieswall/tesis/puntos_fecha_BSHC');

// load and filter the collection
var Landsat8_1 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');
var Landsat8_1_Filtrada =   Landsat8_1.filterBounds(limites_BSHC)
.filterDate('2015-05-4', '2015-08-31')
.filterMetadata('CLOUD_COVER','less_than', 10);

//filter by bounds
var Landsat8_2 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');
var Landsat8_2_Filtrada =   Landsat8_1.filterBounds(limites_BSHC)
.filterDate('2015-08-20', '2015-10-01')
.filterMetadata('CLOUD_COVER','less_than', 10);

//////getting imagery in to variables
//junio
var imagen1 = ee.Image ('LANDSAT/LC08/C01/T1_TOA/LC08_226075_20150807');
var imagen2 = ee.Image ('LANDSAT/LC08/C01/T1_TOA/LC08_226076_20150807');

//septiembre
var imagen3 = ee.Image ('LANDSAT/LC08/C01/T1_TOA/LC08_226075_20150924');
var imagen4 = ee.Image ('LANDSAT/LC08/C01/T1_TOA/LC08_226076_20150924');

//Create mosaics per period
var mosaic_junio = ee.ImageCollection ([ imagen1, imagen2]).mosaic().clip(limites_BSHC);
var mosaic_septiembre = ee.ImageCollection([imagen3,imagen4]).mosaic().clip(limites_BSHC);

//Compute NDVI per mosaic
// Compute the Normalized Difference Vegetation Index (NDVI).
var nir_junio = mosaic_junio.select('B5');
var red_junio = mosaic_junio.select('B4');
var ndvi_junio = nir_junio.subtract(red_junio).divide(nir_junio.add(red_junio )).rename('NDVI_jun');

var nir_septiembre = mosaic_septiembre.select('B5');
var red_septiembre = mosaic_septiembre.select('B4');
var ndvi_septiembre = nir_septiembre.subtract(red_septiembre).divide(nir_septiembre.add(red_septiembre )).rename('NDVI');

print ( ndvi_septiembre);

//Extracting pixel values from sample poligons
var prueba = ndvi_septiembre.gt(poligonos);
//mask pixels
var septiembre = mosaic_septiembre.updateMask(prueba). addBands(ndvi_septiembre);

var ndvif = function(image){
  //add the NDVI band to the image
  var ndvi = ndvi_septiembre.rename('NDVI');
  //isolate the NDVI band
  var quality = ndvi.select('NDVI');
  //get pixels above the POLYGON
  var ndvi01 = quality.gt(poligonos);
  //create a mask from high likelihood pixels
  var ndvimask = image.mask().and(ndvi01);
  //mask those pixels from the image
  return image.updateMask(ndvimask).addBands(ndvi);
};
var collection= ndvi_septiembre.map(ndvif);

It returns an error that says: ndvi_septiembre.map is not a function


texts: OpenLayers GetFeatureInfo for GeoServer Layer Groups

texts:
I've created  a GeoServer layer group. I'm using OpenLayers.Layer.WMS to retrieve the layer group, and this works fine. In addition, I'd like to pull feature info for the layer. For this, I'm using OpenLayer's GetFeatureInfo. But the GetFeatureInfo only returns feature info for only one of the layers in the layer group. Is there a way I can retrieve feature info for both layers in the group? And, is this an optimal way of dealing with layer groups?

A section of my code is shown below:

    var layer_group = new OpenLayers.Layer.WMS(
        "Sample Layer Group",
        "http://localhost:8081/geoserver/wms",
        {layers: "layer_group_name", transparent: true},
        {isBaseLayer: false, visibility: false}
    );
    map.addLayers([layer_group]);

to retrieve the layer group.

And:

    map.events.register('click', map, function (e) {
        var url =  layer_group.getFullRequestString({
            REQUEST: "GetFeatureInfo",
            EXCEPTIONS: "application/vnd.ogc.se_xml",
            BBOX: layer_group.map.getExtent().toBBOX(),
            X: e.xy.x,
            Y: e.xy.y,
            INFO_FORMAT: 'text/html',
            QUERY_LAYERS: layer_group.params.LAYERS,
            WIDTH: layer_group.map.size.w,
            HEIGHT: layer_group.map.size.h});

        window.open(url,
                "getfeatureinfo",
                "location=0,status=0,scrollbars=1,width=600,height=150"
            );
    });

To retrieve the feature info.

My environment setup: GeoServer 2.1.3, PostGIS 1.5, OpenLayers 2.11.

EDIT: The features shown are only for one layer

Increasing FEATURE_COUNT only increases the number of instances(rows) returned. The instances are over 1M. What I'd like to achieve is have one row having feature info for all the layers in the layer group.


texts:
Regarding the GeoJSON representation of a TIN, are there

any specifications or
existing implementations I may not be aware of?

You could probably hack one together using the multipolygon feature type, but this seems inefficient.

Generally speaking, I don't see any surface types (regular or irregular) represented in the GeoJSON world, is there any existing or upcoming work in this area?


texts: Performing Overlay of Two Triangulated Irregular Networks (TIN)

texts:
I am referring to the paper here, section 2.6.1, on the adding and subtracting of two TINs:

  The addition of two TINs can be
  determined exactly and stored into a
  new TIN, because the addition of
  piecewise linear functions again
  yields a piecewise linear function.
  The addition is done by performing an
  overlay of T1 and and T2, there
  are several algorithms for this. After
  this we obtain a subdivision where all
  faces have 3,4,5,6 edges. We now must
  fill in the height information for the
  vertices of the overlay..

Although I can understand every single word out of the passage, I don't know how to carry out the above procedure in practice to obtain the cut/fill of the two TINs.

More specifically, I would like to know how to perform the overlay of two TINs. There are references given at the end of the paper, but I can't have access to them because I'm not inside a university library. So any readily accessible online reference (or code samples) are greatly appreciated!


texts: Why can't I edit polygons once they have been saved to MySQL?

texts:
I have a simple table with an id field and a polygon field in mySql.
Using QGIS I can add new polygons and give them an id. Before saving them I can also edit them.

But once a polygon has been saved it cannot be edited any more: every time I move it, it simply moves back to the original position when I save the edits. Added points disappear when saving.

The same happens with a point field in a different table.

What could be the problem?


texts: How do I get rid of this thin band between raster layers created from TNM?

texts:
After hitting my head against the table for many days, I'm humbly presenting myself to the QGIS community.

I downloaded 4 IMG files (1/3 arc second, elevation data) from The National Map on the USGS website. I then loaded these into QGIS v 2.4 and went through the following steps with each section.

First, I converted the raster from meters to feet using the Raster Calculator (outputting to a GEOTiff file and using the current layer extent).\

Next, I used the DEM(terrain models) tool to create a hillshade raster layer.

I'm assuming my adding of contours doesn't matter.

Why is there a line that forms between the sections? The data is supposed to be seamless between the layers.

Also, my Project CRS is NAD83.


texts: How can I do to export shapefile with symbology.svg COLORED to DXF format?

texts:
I created a shapefile categorized by characteristics, and I assigned to each element a SYMBOL.svg COLORED with Design Effects.
I want to export the shape with the characteristics of the colored symbols.svg that I created and I do:
Project> Export> Export project to DXF ...
but when I open the file.dxf with AutoCAD, it exports only the symbols in black not with the assigned colours.
Do you know how I can do to keep the colours of the Symbology.svg to the shape during the export?


texts: Update attributes with arcpy.updateCursor and updateRow

texts:
I need to update attributes with arcpy.updateCursor and updateRow based on a list of service numbers, but I can't even get one row to update. This is my beginner code.

df= r"GAS FEATURES\SERVICE GROUP\Gas Service Point"
rows = arcpy.UpdateCursor(df)
for row in rows:
   ----if row.getValue("SERVICENUMBER")= '02753-725':
   ----row.SetValue("SubtypeCD",'Multiple Meter Service')
   ----rows.updateRow(row)


texts: How to NOT show decimal values in TileMill when using gdal_contour generated contours?

texts:
Following problem has troubled me for two days now and I can't find a solution:

I am using QGIS to generate contour lines from SRTM DEMs. Export these as shapefiles and load them into TileMill on top of a GeoTIFF.

Problem: I can't seem to get rid of the ".0" end of each height (see picture).

This is the carto.css I use:

#kueste_1000 [zoom>=12]{
  line-width:0.6;
  line-opacity:0.5;
  line-color:#c6b69b;
  ::label{
  text-name: [ELEV];
  text-face-name: @sans_italic;
  text-fill: #c6b69b;
  text-placement: line;
  text-min-distance: 400;
    text-dy:10;
  text-size: 12;
    text-opacity:1;
  }
}

I have different contour shapefiles, generated years ago by dem2topo, labelling these works fine.

I thought the problem might be caused by the attribute label, but changing (in gdal_contour command line) the name of the attribute from ELEV to NAME or anything else doesn't do the trick.

I also tried to change the gdal_contour command line:

gdal_contour ... -i 1000.0

to:

gdal_contour ... -1 1000

but this also didn't work.

The funny part is that QGIS, when setting the flag to label these lines, does so without the .0


texts: dialog for editing QgsGraduatedSymbolRendererV2

texts:
I have a stand-alone pyqgis app, some layers are using the QgsGraduatedSymbolRendererV2 as renderer with customized range and color symbol. I am trying to use an existing Qgs Dialog control to allow users to edit the breakdown and symbol etc of those layers. Out of all of the qgs controls, I found the following as the most probable:

cr = QgsVectorGradientColorRampV2()
ed = QgsVectorGradientColorRampV2Dialog(cr, None)
ed.exec_()

But, I don't think it does what I want as it only allows percent category, not real double values ranges.

So, which editor should I use for QgsGraduatedSymbolRendererV2 type of style? or do I have to make my own?


texts:
This problem is in IE8 and not in Firefox.

My first problem was that IE8 don't support JSON natively and I get the "JSON is undefined". I solved this by including json2 or json3.

Now JSON.parse (both json2 and json3) throws an exception on proj4.js because of a syntax error.

The projection (from the Danish Geodata Agency) used in leaflet is as follow:

 var crs = new L.Proj.CRS.TMS('EPSG:25832',
                              '+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs',
                               [120000, 5900000, 1000000, 6500000],
                               {resolutions: [1638.4, (...), 0.2, 0.1]
 });

I have removed some of the values from {resolutions}.


texts:
I was using OpenDroneMap to generate 2D map and the output comes with corner coordinates like this:

-1.92011886596679688e+02 -1.89522613525390625e+02 1.91508605957031250e+02 1.99496490478515625e+02

This does not look like normal lat/lng data, what format is this?

How should I convert to normal lat/lng?


texts: Are these two Proj4 string treated differently? will any argument get ignored or overrided?

texts:
I have knowledge enough to play around Proj4s string but I have not found any question about this yet.

I have these two Proj4 strings.

+proj=eqdc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=-0.9956,1.9013,0.5215,0.025915,0.009416,0.0011599,-0.00062 +units=m +no_defs

and this one below (+datum=NAD83 is added in front of +ellps=... +towgs84=...)

+proj=eqdc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +datum=NAD83 +ellps=GRS80 +towgs84=-0.9956,1.9013,0.5215,0.025915,0.009416,0.0011599,-0.00062 +units=m +no_deft

In case you wonder, The 7-parameters argument in both strings is from this post.

I wonder if:

In the 2nd string, will +ellps=... +towgs83=... override +datum=NAD83 during transforming projection?
I have two file about the same data projected by 1st string and 2nd string. I am going to transform them to any NAD83-based projection or NAD83 latlong(EPSG-4269). Will Proj4 transform the data projected by 2nd string to the new NAD83-based projection with or without datum shift because I had +datum=NAD83 in the Proj4 string and in front of +ellps=... +towgs84=...?

2.1 After projecting 1st-string-projected data and 2nd-string-projected data to the same NAD83 projection, will or will not I have two same results?


texts: Google Earth Engine: cannot get the properties which are set in the map function

texts:
I tried to set new properties in a map function, but I failed to get these properties outside the map function. The error is:

    ComputedObject (Error)
    Collection.first: Error in map(ID=0):
    Feature.intersection: Parameter 'right' is required.

And the code is:

// double map over the polygons
var mapped = vectors.map(function(feat1){
  feat1 = ee.Feature(feat1);
  var mapped1 = polygons.map(function(feat2){
    feat2 = ee.Feature(feat2);
    var intersection = feat2.intersection(feat1, ee.ErrorMargin(1));
    return ee.Feature(intersection).set({'Intersect': intersection.area().divide(1000 * 1000),'date':'2016-01-01'})
  })
  return mapped1
}).flatten().filter(ee.Filter.gt('Intersect', 0));

var test = mapped.first().get('date');
print(test);

As I set new properties "Intersect" and "date", but I cannot get any of them when I print it.

The full code can be found at:
https://code.earthengine.google.com/e31d1380d8499df42194fa2a5afa08f1

The test code is written by Rodrigo E. Principe
https://gis.stackexchange.com/users/23470/rodrigo-e-principe


texts: Using elevation point of splines from AutoCAD in ArcMap

texts:
I have splines countour lines with Control Point z in autocad . I want use z point values in hecras. So I export polylines of the DWG toArcMap but in the elevation tab the values seems incorrect and zero values.
Normally I have polylines counter lines with elevation so I use this method and it works. But now I have splines countour lines instead of polylines. I tried to convert all splines to polylines but it seems doesn t work.
My aim is using this elevaiton data in order to make analysis in Hecras.
How can I see the correct elevation points in ArcMap?
My AutoCAD file consisting countour lines:
https://easyupload.io/ibi0oy


texts: Shapefile To svg to dxf and back to shapefile

texts:
I need to edit a complex vector polygon layer based on a raster layer. It will be a lot easier for me to do the drawings in an illustration software like Inkscape and then to import it back to QGIS. So I'm planning to export the layers (both the raster and the vector) as svg layers, then to do the edits in Inkscape, and then to save the file as an dxf in Inkscape and finally to import back the dxf to QGIS and converting it back to a shapefile (or geopackage).
The problem obviously is when I import back the dxf to QGIS it's not georeferenced. And if I save the vector layer in QGIS as dxf I can't open it in Inkscape (I don't know why).
Any solutions?


texts: How to save a map service to a local drive.. Forest Service Ranger Districts

texts:
I am using a map service in a project. I would like to find out if it is possible to save it to my local computer.

This is what I am looking at:

Server: http://apps.fs.fed.us/ArcX/services
Name: EDW/EDW_RangerDistricts_01
Map Service Type: Not Cached


texts: How to iterate folder in a model that has a sub-model (iterate raster) in ModelBuilder

texts:
I have multiple folders that contain a number of raster files. What I would like to do is iterate through every folder and resample each every raster file from 1km to 30m.

Since only one iterator is allowed in one model, I figured importing the submodel (iterate raster) into another model (iterate workspace).
In the 'Resample30m' model, I created a folder variable and defined it as a model parameter.
I chose 'folder' for Workspace type in the workspace iterator, and tried to connect the iterate workspace's output (Workspace) with Resample30m model, it was impossible to do it. Only the environments and precondition option popped out.

I'm guessing the iterate workspace doesn't access the inner composition and that's the one raising the problem.

Is there any way that I can achieve what I want in ModelBuilder?


texts: Using LineFillSymbol in polygon layer from ArcObjects?

texts:
I'm creating a symbol hatch:

var lineFill = new LineFillSymbolClass
{
    Angle = 45,
    Offset = 0,
    Separation = 0.1,
    Outline = new SimpleLineSymbolClass { Style = esriSimpleLineStyle.esriSLSSolid, Color = Color.FromArgb(215, 15, 245).ToRgbColorClass(), Width = 0.5 },
    LineSymbol = new SimpleLineSymbolClass { Style = esriSimpleLineStyle.esriSLSSolid, Color = Color.FromArgb(215, 152, 245).ToRgbColorClass(), Width = 0.5 }
};
((IGeoFeatureLayer)featureLayer).Renderer = new SimpleRendererClass { Symbol = lineFill };

Instead of hatch, makes the background simple fill.
What am I doing wrong?


texts: Raster reprojection from EPGS:4269 to EPSG:4326

texts:
I downloaded dem rasters for USA and they use EPSG 4269 CRS. Now, I want to sample raster values (elevation) at some points that use EPGS:4326 CRS.
Do I need to reproject rasters to 4326 before sampling values?


texts: Removing points which fall outside of selected polygons using QGIS?

texts:
I have a project of finding the safest and most convenient place to cycle in London. There is a polygon layer(the yellow polygons) of boroughs that have lower collision rate; and a points layer showing all the parking space (green points).

I wanted to show only points in yellow polygon, I have tried to use vector--> Geoprocessing tools --> clip, I used points layer as the input layer, and yellow polygon layer as the clip layer. But the result removed all the points on map.

Can someone tell which step did I do wrong ?


texts: Geowebcache disk quota is different from tile layers page

texts:
This is my disk quota page and following is my tile layers list page:

And I could not found my tiled images in server data directory.

I think my tiles are not saving on disk. What is the problem that you think?


texts: Average of raster elevation based on polygon boundaries

texts:
I am attempting to extract and calculate the average height of buildings using LiDAR converted to raster elevation data. I've created a DSM to cover the entire area and I am trying to extract just the buildings out, I need the average height based upon the DSM elevation pixels. How do I calculate this?

I've attempted clipping these images out according to the polygons (which created thousands of individual rasters) but still I am unsure how to get the average within the bounds of the footprint.


texts: Setting up GeoServer Image Mosaic plug in to pick up images from directory dynamically?

texts:
I have set up Image Mosaic to load GeoTIFF images from a directory and it is all good but when I add more GeoTIFF files to the directory they do not show up in the mosaic unless I recreate it.

Is there a way to configure Image Mosaic to pick them up automatically?


texts:
I'm trying to draw my geometery on worldwind but it keeps saying

  WARNING: Unspecified coordinate system

and draws a default blue line in the middle of nowhere.

The procedure is that I populate my db with a geometery like this:

INSERT INTO newroutesegments (name, geom) VALUES ('RouteSegment1', ST_GeomFromGeoJSON('{ "type": "MultiLineString", "coordinates": [ [ [100.0, 0.0], [101.0, 1.0] ], [ [102.0, 2.0], [103.0, 3.0] ] ], "crs":{"type":"name","properties":{"name":"EPSG:4326"}} }'));

At this point I checked my db and the geom is there and the SRID is 4326 which is good! However when I query this and want to draw it on world wind it gives me the same error. This is my geometery in case you want to see it:

0105000020E610000002000000010200000002000000000000000000594000000000000000000000000000405940000000000000F03F010200000002000000000000000080594000000000000000400000000000C059400000000000000840


texts: OpenLayers: How to re-add a feature that draws line segments of a tracked route?

texts:
My app records tracks (e.g. a bike route). There is a start and a stop button, which start and stop the recording of such a track. After a recording has been stopped, the user can start a new recording again. And here start the troubles: When recording again, the LineString ist not drawn again.
Here is the relevant code:
When the recording is started, these things happen;
// in Vue methods
// "this.trackFeature", "this.trackLayer" etc are defined in Vue's "data"
// and setup in the "mounted()" lifecycle

startRecording () {
  // Create / add trackLayer if there isn't any
  if (!this.trackLayerAdded) {
    this.trackFeature = new Feature({
      geometry: new LineString([])
    })

    this.trackSource.addFeature(this.trackFeature)
    this.trackLayer.setSource(this.trackSource)
    this.trackLayer.setStyle(this.trackStyle)

    // Check if "this.trackLayer" has already been added, otherwise
    // openlayers' "reassign" error (58) is thrown
    let trackLayerExists = null
    const layerCollection = this.olMap.getLayers()
    layerCollection.forEach(l => {
      if (l.ol_uid === this.trackLayer.ol_uid) {
        trackLayerExists = l
      }
    })

    // If the "trackLayer" was not found, add it
    if (!trackLayerExists) {
      this.olMap.addLayer(this.trackLayer)
    }

    this.trackLayerAdded = true
  }
}

Since the line style is defined with { radius: 30, fill: [_yellow in rgba_] }, a big yellow point appears on the screen on my desktop computer (where the geolocation` does not move).
When the recording is stopped and saved, this happens:
// in Vue methods
deleteTrack () {
  this.trackSource.removeFeature(this.trackFeature)
  this.trackSource.clear()
  this.trackFeature = null // needed?

  // // Do I have to remove the layer, too?
  // this.olMap.removeLayer(this.trackLayer)

  this.trackLayerAdded = false
},

The yellow line segment disappears.
Now, when I click on "start" again to record a new route, the yellow LineString does NOT appear anymore, which is, what I would have expected.
I cannot find my mistake. What am I doing wrong? What do I have to do that each time I click on "start record" a new LineString will be recorded?


texts: How to calculate distance from POINT to LINESTRING in R using `sf` library and get all POINT features from LINESTRING where distances were calculated?

texts:
I'm trying to calculate distances from a set of POINT features (1, 2, 3, 4, 5, 6 & 7) to a LINESTRING feature using sf library in R. I could calculate all the distances using the function st_distance but I can't get which POINT feature from LINESTRING feature were used by the function to calculate those distances:

First example: if the minimum distance between POINT "2" and the LINESTRING is the distance between POINT "2" and the POINT "A" in the LINESTRING, how can I get POINT "A"?
Second example: if the minimum distance between POINT "4" and the LINESTRING is the distance between POINT "4" and the POINT "X" that is located between POINT "B" and POINT "C" in the LINESTRING, how can I get POINT "X"?

Note: I know it's possible to achieve this using geosphere::dist2Line function but at the moment of this question geosphere package don't have compatibility with simple features objects. So I'm trying to solve this issue using sf without using Spatial objects from sp.

Below is a reproducible code example and a figure of the problem:

# Load Libraries ----------------------------------------------------------

library('sf')

# Test data ---------------------------------------------------------------

points.df <- data.frame(
  'x' = c(-53.50000, -54.15489, -54.48560, -52.00000, -52.57810, -49.22097, -48.00000),
  'y' = c(-38.54859, -39.00000, -38.80000, -38.49485, -38.00000, -38.50000, -37.74859)
  )

line.df <- data.frame(
  'x' = c(-52.53557, -52.00000, -50.00000, -48.00000, -46.40190),
  'y' = c(-41.00000, -40.60742, -40.08149, -40.82503, -39.00000)
  )

# Create 'sf' objects -----------------------------------------------------

points.sf <- st_as_sf(points.df, coords = c("x", "y"))

st_crs(points.sf) <- st_crs(4326) # assign crs
points.sf <- st_transform(points.sf, crs = 32721) # transform

line.sf <- st_sf(id = 'L1', st_sfc(st_linestring(as.matrix(line.df), dim = "XY")))
st_crs(line.sf) <- st_crs(4326) # assign crs
line.sf <- st_transform(line.sf, crs = 32721) # transform

# Plots -------------------------------------------------------------------

xmin <- min(st_bbox(points.sf)[1], st_bbox(line.sf)[1])
ymin <- min(st_bbox(points.sf)[2], st_bbox(line.sf)[2])
xmax <- max(st_bbox(points.sf)[3], st_bbox(line.sf)[3])
ymax <- max(st_bbox(points.sf)[4], st_bbox(line.sf)[4])

plot(points.sf, pch = 19, col = "#53A8BD", xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin,xmax), ylim = c(ymin,ymax), graticule = st_crs(4326), axes = TRUE)
plot(line.sf, col = "#C72259", add = TRUE)
text(st_coordinates(points.sf), as.character(1:7), pos = 3)
text(st_coordinates(line.sf), LETTERS[1:5], pos = 1)

# Distances ---------------------------------------------------------------

distances <- st_distance(x = points.sf, y = line.sf)
print(distances)

Units: m
         [,1]
[1,] 262727.8
[2,] 256710.2
[3,] 292476.9
[4,] 223153.4
[5,] 291143.6
[6,] 188868.4
[7,] 198670.4


texts: How to write unique ID(s) of polygons in layer A intersecting polygons in layer B to field in layer B?

texts:
I'm working in ArcGIS 10.1. Suppose that I have two polygon layers, A and B.

Polygons in layer B can have one of a few relationships to those in layer A:

They may be completely contained by polygons in layer A
They may intersect (be shared by) multiple polygons in layer A
They may be completely contained but coextensive with those in layer A

I know of a number of ways to visually identify, or select for, these relationships, but I'm wondering if anyone has some idea of how I could write, to the attribute table of B, values from the unique ID of layer A based on these relationships?

FWIW, I'm working with the relationships among counties, county subdivisions, and places in the TIGER hierarchy -- i.e., where A = County Subdivisions and B = Places, and in turn where A = Counties and B = County Subdivisions.  It's possible that I'm neglecting some of the potential relationships between A and B in the way I've laid out the list above.

Thanks in advance for any help you can provide!


texts:
When using the refresh() function, is it possible to only refresh part of the map canvas, instead of refreshing all of it?

I am using C++ but can use any code snippet even if it is in Python.


texts: Problem with query in a loop in QGIS Script

texts:
I want to run the following script in QGIS

from qgis.core import *
from qgis.utils import *
from PyQt4.QtCore import *

inLayer = iface.activeLayer()

with open('iso_codes.txt') as f:
  codes = f.readlines()

for isoCode in codes:
  #exp = QgsExpression('"ISO_CODE" = \'AF\'') # THIS WORKS
  exp = QgsExpression('"ISO_CODE" = \'%s\'' % isoCode) # THIS DOES NOT WORK
  request = QgsFeatureRequest( exp )
  fipsCode = ''
  for feature in inLayer.getFeatures( request ):
    fips = feature['FIPS_CODE']
    if fipsCode != fips:
      fipsCode = fips
      print '%s, %s' % (isoCode, fipsCode)

Using '"ISO_CODE" = \'AF\'' as argument the script works till the end.

Using '"ISO_CODE" = \'%s\'' % isoCode as argument the script doesn't go into the for loop.

I've tried to use inLayer.setSelectedFeatures([]) at the end of the for loop, but this did not help.


texts:
I am a proficient GIS and Python user. I would like to ask if any open source tools exist which allow for rendering of geospatial data on a website. I would prefer to avoid tools which require in-depth understanding of underlying client/server architectures or learning new programming languages such as JavaScript or HTML.

Do any such tools exist? Ideally I would like to be able to rapidly deploy a web viewer containing a map and some layers.


texts: Adjust a complete shapefile of multiple points on ArcMap taking other few as reference?

texts:
I received a shapefile with georeferenced trees, but they didn't use a GPS to get the coordinates (not sure how they did it, some traditional stuff maybe?). So, the position doesn't match exactly on the field. I need to correct this for planning some plots. I mapped myself some of the trees on the field and compared to the first shapefile, they seem to have a semi-consistent displacement at NE and the distance varies between 70-80 m.

Is there some way to correct this? I can afford accuracy to 10 m

If so, how many trees would be needed to map on the field for reference? The first shapefile has approx. 3000 points


texts:
I am working with the Flow Accumulation tool in ArcMap 10's spatial analyst hydrology tools. This tool calculates the total number of raster cells that flow into a given cell. However, for my purposes, I would like to restrict how far uphill the accumulation originates.

IE: I would like to see how much flow there is into a cell within 1 km instead of all the way to the ridge-line.


texts: Identify layers in an .mxd that are network elements feature types

texts:
I need to create new feature classes from all layers in an .mxd, and also create a new .lyr file for each of the newly created feature classes that retains the symbology of the original map layer. The script works for most files, but fails on layers that have network element feature types- Junctions and Edges.

http://desktop.arcgis.com/en/arcmap/10.3/manage-data/geometric-networks/essential-geometric-networks-vocabulary.htm

It creates the new feature class on these layer types, but fails on creating the layer file. This is because the original map layer is of feature type: "simple junction" for instance, and the newly created feature class doesn't retain that type. It is now feature type: "simple." When I try to apply symbology I then get the error- 000968: The symbol layer does not match the input layer.

I either need to:

Find a way to maintain the feature type if it's a network element layer. I've tried using the FeatureClassToFeatureClass tool, and CopyFeatures tool, and neither has worked.

or

Identify the layers that have a network element feature type so that I can deal with them in another way. I have tried using both the describe and layer properties, and haven't found a way to identify these layers.

Relevant portion of code causing the error:

if not arcpy.Exists(fc_outname):
    arcpy.AddMessage("Creating feature class: {}".format(fc_outname))
    out_path, out_name = os.path.split(fc_outname)
    arcpy.FeatureClassToFeatureClass(clipped, out_path, out_name)

#Create .lyr file
if not arcpy.Exists(lyr_outname):
    arcpy.AddMessage("Creating .lyr file: {}\n".format(lyr_outname))
    lyr = arcpy.MakeFeatureLayer_management(fc_outname, lyr_disp)
    arcpy.ApplySymbologyFromLayer_management(lyr, v)
    arcpy.SaveToLayerFile_management(lyr, lyr_outname)

Using ArcGIS 10.5


texts:
I have made a buffer on the area shown. If you closely in zoom there is a line beneath the buffer. I am using ST_Within for the comparison. There are two linestrings in my area. I am making buffer of about 10 meters on one area and comparing it with the other linestring using ST_WITHIN. The linestring is divided into different segments. ST_Within checks that two geometries should be completely inside. If you can see closely in the picture I only want the linestring which is inside the buffer but ST_within ignores the complete segment because small portion of line is outside buffer and Moreover, I also dont want this extra line which is outside buffer. Is there a way to compare these lines so that I can get only that segment which is inside buffer?. I cannot use ST_Intersects because then I will have other extra intersecting lines as well which I dont want.
Here is my query:

Create Table CombinedTable As select area1.Geometry AS Geometry from Area1 area1,Area2 area2 where
St_Within(area2.geom,ST_Buffer(SetSrid(area1.geometry,3857),10))


texts: Parts needed to run Geoserver/PostGIS from a program like Microsoft Access

texts:
I have a fair bit of experience in GIS and Python, but am brand new to developing on the web.

I have rented some space on a server and have a PostgreSQL/PostGIS database as well as Geoserver on it.  The other side of the coin is a Microsoft Access database (which will later be a web based application, but that's later) that has a bunch of data related to the Spatial data (polygon layers) up on the server.

What components do I need to develop so that the Access database can request specific data from Geoserver and have it display on a webpage (that is seen by the user through the Access program).  I also want the Access program to be able to request a KML layer that will display in Google Earth.

For the first option I think that I need to use the OpenLayers GetMap examples and build a map from that.  The map will need to be able to get all of a particular layer, as well as one polygon from the same layer displayed differently so that it stands out.  I remember that there is a way to integrate Google Maps as a base layer.  Are there other/better options?

For the second option (use the Access database to query Geoserver and display the result in Google Earth) I'm at a loss.  It seems like it can be done, but how?

I just need to be pointed in the right direction as most of these technologies are new to me and I'm a bit lost.

Thanks a lot!


texts: How do I create a buffer distance from the edge of the polygon in qgis? The buffer already created calculates the distance form the center of the polygon. For my analysis, I'll need a hint for calculating a buffer form the edge of the polygon.


texts: ArcPy FeatureClassToFeatureClass_conversion invalid SQL statement

texts:
I am trying to iterate through a list of where clauses as such:
for vc in vcList:
    print (vc)
""""D_Zone" =  1"""
""""D_Zone" =  2"""
""""D_Zone" =  3"""
""""D_Zone" =  4"""
""""D_Zone" =  5"""

for each iteration I want to use a for loop to execute aFeatureClassToFeatureClass_conversion using the iterable as the where clause:
for vc in vcList:
    oName = base + "_" + str(count)
    arcpy.conversion.FeatureClassToFeatureClass('Origins_Geocoded', oPath, oName, vc)
    count += 1

but I am getting an invalid SQL error:
An invalid SQL statement was used. [SELECT * FROM Origins_Geocoded WHERE """"D_Zone" =  1"""]

The value being queried is short numeric. Seems to be consistent with the first example from Esri:
https://pro.arcgis.com/en/pro-app/arcpy/get-started/specifying-a-query.htm
What must I do to construct a proper where clause?


texts: Change attribute table name in Tiger file

texts:
I'm fairly new at QGIS, though pretty experienced with ArcMap. So I thought I would learn it with a fun data set, a Census bureau tiger file. Specifically the Public Use Microdata Area (PUMA) which can be found at.

https://www.census.gov/geo/maps-data/data/tiger-data.html

Each of the attributes has a alphanumeric name which corresponds to the metadata with the actual descriptions. When I try to open the metadata QGIS slows down until it's unusable. How can I transfer or best leverage the metadata to know what the attributes are. I know it's a broad and somewhat simple question, but I really cannot find a good way to do this.


texts: Storing results of an iteration Google Earth Engine

texts:
I am trying to get summary statistics of pixels within a buffer region. For this, I am selecting a geo-coordinate that I am interested in, defining a 10 km radius around it and counting the number of pixels which are no data, not water and water. This is the only information stored in the monthly water history data.
I have a question regarding the following code. I don't understand how I can store the results I get from each iteration over the separate months in the image collection named water. Any suggestions on how to store results?

// 1. Import water data and filter months
var gsw = ee.ImageCollection("JRC/GSW1_1/MonthlyHistory");
var water=gsw.filterDate('2018-01-01','2018-12-31');
// 2. Sample coordinate
var test=ee.FeatureCollection(ee.Geometry.Point([38.8,7.97]));
//3. Define buffer region for the coordinate: 10,000 kms
var buffer=test.map(function(feat){
  return(feat.buffer(10000));
});
// 4. Create an empty list to store the count of pixels of no data, not water and water in each month for the same buffer region
var month=ee.List([]);
var extract=water.iterate(function(image,first){
    var image0=image.updateMask(image.select('water').eq(0));
    var image1=image.updateMask(image.select('water').eq(1));
    var image2=image.updateMask(image.select('water').eq(2));

    //5. create the count of pixels of the three categories 'nodata', 'notwater','water'.
    var reduce0=image0.reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: buffer.geometry(),
      scale: 30,
    });
    var reduce1=image1.reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: buffer.geometry(),
      scale: 30,
    });
    var reduce2=image2.reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: buffer.geometry(),
      scale: 30,
    });
    var reduce=image.reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: buffer.geometry(),
      scale: 30,
    });

    var nodata=reduce0.get('water');
    var notwater=reduce1.get('water');
    var water=reduce2.get('water');
    var total=reduce.get('water');

    return month.add(nodata, notwater,water);
  });


texts: Creating Stored Procedure in versioned SDE_SQLServer geodatabase using DDL?

texts:
I've almost used SQL server_SDE geodatabases for 3 years in GIS projects and I don't know how can I create a stored procedure in geodatabase.I use ESRI case tool to create a geodatabase in ArcCatalog(using xml file as described in ESRI documentations),so I never created a geodatabase with DDL scripts,But i would like to know (and do) creating a geodatabase using DDL commands in SQL Server(while using SDE and versioning) and adding stored procedure to this geodatabase.As I know,sde creates some metadata tables and breaks data of versioned tables into some other tables(like a_xx and f_xx,for example a89,f89) so that  some records of table is saved in these helping tables,so that I don't know how can I find the name of these helping tables and basically do i have to work with them in stored procedure or not.

Is there any document or guide line to show how can create a geodatabase using DDL command in an environment like SQL server management studio?and how can create stored procedures in geodatabase and how can call it?


texts:
I successfully created a shapefile from OpenStreetMap with the OpenStreetMap tool in QGIS 2.2. When I try to create a polyline layer, it says that there are 8000 values in the db file for the property highway, but after the import I have only 6000 in the attribute list.

The problem is that I don't have all streets available for a given city and because of that I am unable to draw parcels where the streets are their borders. Is that problem normal or am I missing something?


texts: How can I buffer a mulipolygon only on the coastline?

texts:
How can I achieve something similar to the below image? (Apologies for the crudely drawn diagram) Either in QGIS or PostGIS.

I have a series of station grounds and I want to retain the inner boundaries, but I want the coastline to buffer out 1 mile (distance for example, could be any distance). Also I want to split an estuary equidistant between polygons.


texts: Changing census income data from string to number in QGIS

texts:
I am trying to use Census Block group income data that I retrieved from American Fact finder. When I join the income data to the block group shapefile in QGIS it imports as a string.

How can I convert this in to an integer?


texts: Generate a matrix of change from shoreline vector time-series data

texts:
I have shoreline vector data that corresponds to property ownership, and time series data referring to attributes such as erosion levels, shoreline class (ornamental, natural, regenerative), etc.

At this point I have two years of data, and I'm looking to generate a matrix of change, that describes from-to change, and no change (Natural shoreline to Ornamental, or Natural to Ornamental, as examples).  With 4 shoreline classes, there are 16 possible change scenarios.

Any suggestions, I primarily use Arc10, but I am also enjoy open source solutions.

In advance, thanks for your time!


texts:
I'm trying to use arcpy.Statistics_analysis using a feature layer as the input and a long int field ("MaxOID5" for the field summary so I can get the max value.

arcpy.env.workspace = r"C:\data\mydata.gdb"

fc = "airplanes"

arcpy.MakeFeatureLayer_management(fc, "fc_layer")

arcpy.Statistics_analysis("fc_layer", "Max_Table", "MaxOID5", "MAX")

When I run the script I get the following error:

ERROR 000800: The value is not a member of SUM | MEAN | MIN | MAX | RANGE | STD | COUNT | FIRST | LAST.

The help on the error says that "The keyword used is not a member of the list of acceptable keywords" but I'm choosing it directly from the auto generate parameters in the arcmap python command window.


texts: QGIS 3.10 open data source manager: How to load XYZ tiles

texts:
I have a thunderforest background map for my project for which I had not given an API key.  I now want to publish some maps and need to add the key to remove the overlay.
The problem is that I can't figure out how to change the urls for the layers.  I think I loaded it originally from the data manager => XYZ Tile menus (this matches what I found searching) but in the current version of QGIS that option is not listed.

I assume I could hack the project file and add the key to the url there but I would rather do this properly!


texts: Would QGIS 2.8 print composer_templates work in QGIS 2.14?

texts:
I have created some maps using print composer of Qgis version LTR 2.8.9. and saved them as templates.

I am planning to upgrade to 2.14.1 LTR version and I will completely uninstall Qgis in order to have a proper clean install of the new version. I will backup the folder composer_templates from c:\users\user\.qgis folder. Will the templates be ok with the new version? Has any one tried the same even with other versions?


texts: Automatic conversion to surface reflectance (SCP)

texts:
I am practicing for SCP part by using Qgis 2.10.1.
I already download desired LANDSAT image from SCP>Download Landsat.

However, in Pre-processing part, after click Select directory button, the list of bands do not appeared in Metadata part for me to do the Perform Conversion.

These are the data download from SCP>Download Landsat

Is there something wrong with my data perhaps?


texts: Create .lyr file from open .mxd layer without removing layer

texts:
I have a tool that is meant to clip and reproject every layer in an open .mxd, and then create a new feature class and .lyr file with original symbology from the result, not adding either to the .mxd. I would like the new .lyr file to have the same display name as the original, but when I write it this way it removes the original .lyr from the map document, and also fails on importing symbology.

I'm using ArcGIS 10.5.

Is is possible to create a new .lyr with the same name as the open .mxd .lyr, then import symbology, and not have it crash remove the original .mxd .lyr?

My current workaround is to append "_clip" to the arcpy.MakeFeatureLayer_management out_layer parameter as follows:

if not arcpy.Exists(lyr_outname):
    if not os.path.exists(out_lyr_fldr):
        os.makedirs(out_lyr_fldr)
    arcpy.AddMessage("Creating .lyr file: {}\n".format(lyr_outname))
    lyr = arcpy.MakeFeatureLayer_management(fc_outname, i.name + "_clip")
    arcpy.ApplySymbologyFromLayer_management(lyr, i)
    arcpy.SaveToLayerFile_management(lyr, lyr_outname)

This works, but people don't want the "_clip" in the display name.

If I try it using the original name:

if not arcpy.Exists(lyr_outname):
    if not os.path.exists(out_lyr_fldr):
        os.makedirs(out_lyr_fldr)
    arcpy.AddMessage("Creating .lyr file: {}\n".format(lyr_outname))
    lyr = arcpy.MakeFeatureLayer_management(fc_outname, i.name)
    arcpy.ApplySymbologyFromLayer_management(lyr, i)
    arcpy.SaveToLayerFile_management(lyr, lyr_outname)

it first removes the layer its working on from the .mxd, and then crashes on import symbology.

Failed to execute. Parameters are not valid.
ERROR 000968: The symbol layer does not match the input layer
Failed to execute (ApplySymbologyFromLayer).

If there is no direct way to achieve this, I thought maybe I can just search the output folder and change everything after it has been run. I can see the display name using describe nameString:

workspace = r"C:\..."
lyrs = []

walk = arcpy.da.Walk(workspace)

for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        lyrs.append(os.path.join(dirpath, filename))

for lyr in lyrs:
    desc = arcpy.Describe(lyr)

    if "_clip" in desc.nameString[-5:]:
        print desc.nameString

But is there a way to change the nameString?


texts: Altering composer label items in QGIS with Python

texts:
I'm only new to QGIS and Python so this may exists already, and hopefully be easy to do.

What I'm hoping to do sometime in the future is build a map book generator using QGIS as the base. In order to do this I need to be able to select all the labels in the composer and replace some tags, so something like $NAME$ would replace with a variable the user has set at runtime.

I have been having a look at the QGIS APIs [and source] but just can't seem to get it to do what I need. The code I have so far:

composerList = self.iface.activeComposers()
if(len(composerList) < 1):
      return

composerView = composerList[0]
composition = composerView.composition()
if(composition is None):
      return

for item in composition.selectedComposerItems():
//Need to check for QgsComposerLabel type and update the text.

I'm stuck when it comes to checking the type of the current item in Python, selectedComposerItems() just returns each item as a QgsComposerItem and calling displayText() fails because it doesn't exists on the type QgsComposerItem.

I have tried the following code but it crashes QGIS when it hits a legend or scale bar:

for item in composition.selectedComposerItems():
    try:
        item.__class__ = qgis.core.QgsComposerLabel
        if item.displayText() == "$NAME$":
            print "I found one"
            item.setText("Hello World")
    except:
        print "Fail"

The idea was to just try setting __class__ on every item as qgis.core.QgsComposerLabel and if it throws an exception then I know it's not a label, but at the moment if it hits a legend or scale bar it crashes QGIS which is not good.


texts: ERDAS IMAGINE: Spatial Modeler Volume Change

texts:
Goal: sum of volume change of two DEM; time 1 - time 2

How to insert an AOI/query for certain area

How to get visual change of loss and gain?


texts:
I have an aermod output data which has angle (degrees), distance (m) and pollutant concentration from a point (x, y).

How can I create points based on the angle and distance with concentration as attributes in ArcGIS?


texts:
I have a huge problem that I can't solve even after trying everything I can think of. I want to merge 4 rasters as you can see in this image
I have Landsat 8 images, to not have a problem with bands, I have just used Band4 of Landsat 8 and I didn't change their projection EPSG4326 (the same projection for the 4 areas) and for NO Data I have used 0;

and after I have done the operation of merging them, I have found a problem; only 2 rasters are merged and the two have changed their position like this;


texts: PyQGIS: "for" loop with nested "if ..elif ..." condition

texts:
I want to update an existing field according to another, it means the equivalent of a CASE WHEN .... THEN ..... ELSE .... END in SQL , but I want to know how to write this same update in Python for QGIS:

I wrote the following code:

#TEST : mettre à jour un champ existant avec boucle
from PyQt4.QtCore import *
from qgis.core import *
vl = iface.activeLayer()
vl.startEditing()

#step 1
vl.updateFields()
idx = vl.dataProvider().fieldNameIndex( 'ETIQUETTE' )
adduc = vl.dataProvider().fieldNameIndex( 'MODE_ADDUC' )

#step 2
for f in vl.getFeatures():
    if f[adduc] == 'AERIEN':
        f[idx] = 'CIEL'
    elif f[adduc] == 'FACADE':
        f[idx] = 'MUR'
    elif f[adduc] == 'SOUTERRAIN':
        f[idx] = 'TERRE'
    else:
        f[idx] = 'RIEN'
    vl.updateFeature( f )

vl.commitChanges()

but but I got a Python error :

Traceback (most recent call last):
  File "<input>", line 2, in <module>
KeyError: '-1'

Can you help me, please ?


texts: Snapping multiple nodes of a polygon to the nodes of a line?

texts:
I am working with vector polygons that have derived from hand drawn areas from when mapping was carried out with tracing paper over paper maps. Therefore many of the polygons don't exactly line up to the boundaries on the ground.

Another dataset of vector lines I have is accurate and so I have been manually snapping nodes of the polygons to the nodes of the lines. But this is a lengthy process.

Here is a screenshot...

The green line is the accurate dataset. You can see how the polygon boundaries (orange) are roughly correct. Bottom left I have started to edit the nodes, with the snapping options set up, and this works to manually move the nodes of the polygons to snap to the nodes of the line. (Yes, I could have enabled topographical editing to move both the polygon boundaries together, and I do do that).

Is there a method for snapping multiple nodes of the polygon to the nodes of the lines, or is there perhaps a better method of doing this? Using Python perhaps? I am looking for something akin to the 'Enable tracing' function (which is great) but for moving polygon boundaries.

I am aware that this may in places be a duplicate as there are many similar questions, but I haven't quite found anything that answers my particular situation? In particular I have looked at...

QGIS how to align existing polygon borders

How to snap polygons to lines?

How to snap lines to nearest point?

https://gis.stackexchange.com/search?q=how+to+snap+many+polygons+to+lines


texts:
It seems an easy question, but I can't find the Python equivalent of R's Dissaggregate:

  Disaggregate a RasterLayer to create a new RasterLayer with a higher
  resolution (smaller cells)

I expected it in GDAL but no luck so far.


texts: How to draw a custom map on top of base layer

texts:
I'm looking for the easiest way to draw a custom layer and export it as a tile cache.   Basically what I would like to be able to do is load a base layer, and then draw objects and annotations on top of it and export these objects and annotations as a raster (GeoTIFF) than I can them run through MapTiler.  I'd like my custom map to include both polygons, circles, and text, sort of like using a georeferenced Adobe Illustrator.   What's the best process to do this?

I'm thinking QGis since I have experience there, but maybe there are other tools where this is more straightforward


texts: Get the GEBCO 2014 box from latitude and longitute

texts:
I am trying to get the information from

http://coastal-mapping.eu/

It is a great viewer that shows the depth or altitude of a coordinate and the contour if you select tools and draw a line.

The issue is that it uses a grid system instead of latitude longitude: for instance calling to 28.4066N, 16.3033W which creates the following call:

http://wxs-coastal-mapping-eu.aw.atos.net/wms/r/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=BODC-GEBCO&LAYERS=BODC-GEBCO&TILED=false&INFO_FORMAT=application%2Fjson&FEATURE_COUNT=1&I=29&J=0&WIDTH=256&HEIGHT=256&CRS=EPSG%3A3035&STYLES=&BBOX=1039062.5%2C1741699.21875%2C1042480.46875%2C1745117.1875' -H 'Accept: */*' -H 'Accept-Encoding: gzip, deflate' -H 'Accept-Language: es-ES,es;q=0.8,en-US;q=0.5,en;q=0.3' -H 'Connection: keep-alive' -H 'Host: wxs-coastal-mapping-eu.aw.atos.net' -H 'Origin: http://coastal-mapping.eu' -H 'Referer: http://coastal-mapping.eu/' -H 'User-Agent: Mozilla/5.0 (Windows NT 10.0; WOW64; rv:47.0) Gecko/20100101 Firefox/47.0

Where creates an special BBOX based on the GEBCO system.

BBOX=1039062.5,1741699.21875,1042480.46875,1745117.1875

Is there a way to convert latitude, longitude into that BBOX system?


texts: Export Sentinal 2 TIFF image from SNAP with all the bands intact?

texts: I have some Sentinal 2 data from which I need an image of a small area. How do I export that particular image with all the 13 bands intact to a TIFF image using SNAP? I can download 3 bands by viewing RGB and then exporting but I need all the bands intact.


texts: Not able to find some of the gdal drivers when using the java bindings as opposed to gdal by itself

texts:
I am new to the gdal scene so there might be a simple explanation for my problem.  When executing ogrinfo from the command line there is a large list of supported formats.

    C:\OSGeo4W\bin>ogrinfo.exe --formats
    Supported Formats:
    "FileGDB" (read/write)
    "ESRI Shapefile" (read/write)
    "MapInfo File" (read/write)
    ...

However, when I run the equivalent command through the java bindings I get the following output:

    C:\OSGeo4W\bin>java -classpath ..\lib\gdal.jar;. ogrinfo --formats
    Supported Formats:
    "ESRI Shapefile" (read/write)
    "MapInfo File" (read/write)
    ...

Notice that the "FileGDB" driver doesn't appear when accessing gdal through the java bindings.  Why are some of the gdal drivers not available with the java bindings? Is it a problem with my classpath and some of the driver dlls locations?


texts: QGIS composer - How can I eliminate striping across transparent raster layers?

texts:
I have downloaded QGIS dev version 1.7.0-114 and I have come across some weird behaviour in the composer.

I create a map with a hillshade and a transparent DEM ontop of it, pretty standard stuff. When I export this map to an image format I get two strips across my output, this happens for a number of formats (png, jpg, bmp, tiff) and the stripes seem to be in the same place everytime. I have attached an example below (sorry for the size of the attachment).

The stripes appear in the same position on the page every time and this happens with different raster datasets, so I have ruled out the data.

Also depending on the page size you may get 1 stripe (A5 as shown above), 2 stripes (A4) or no strips (custom dimensions h100 x w100).

Is this something that can be overcome via some settings which I may have overlooked? Or might this be something funny in the software?


texts: Mapping over list of dates in Google Earth Engine

texts:
My function fnc returns as Image but on using the ee.ImageCollection(dates.map(fnc)), it returns an error ImageCollection (Error)
ImageCollection.fromImages: Attempt to create an ImageCollection with non-image elements.
 in the code below.

var Date_Start = ee.Date('2015-05-01');
var Date_End = ee.Date('2017-12-01');
var Date_window = ee.Number(30);

// Create list of dates for time series
var n_months = Date_End.difference(Date_Start,'month').round();
var dates = ee.List.sequence(0,n_months,1);
var make_datelist = function(n) {
  return Date_Start.advance(n,'month')
}
dates = dates.map(make_datelist);

var fnc = function(d1) {
  var start = ee.Date(d1);
  var end = ee.Date(d1).advance(1,'month');
  var date_range = ee.DateRange(start,end);
  var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filterDate(date_range)
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))

  return(S1.first())
}

var mt = ee.ImageCollection(dates.map(fnc));
print(mt)

How can I get the image collection back from the mapped function?

EDIT
After trying further, I am now getting a similar error while mosaicing two images (this time from Landsat 5) and then returning the mosaic. Code is below:

var p1 = ee.Geometry.Point([103.521,13.028]);
var p2 = ee.Geometry.Point([105.622,13.050]);

var Date_Start = ee.Date('2000-05-01');
var Date_End = ee.Date('2007-12-01');
var Date_window = ee.Number(30);

// Create list of dates for time series
var n_months = Date_End.difference(Date_Start,'month').round();
var dates = ee.List.sequence(0,n_months,1);
var make_datelist = function(n) {
  return Date_Start.advance(n,'month');
};
dates = dates.map(make_datelist);
var fnc = function(d1) {
  var start = ee.Date(d1);
  var end = ee.Date(d1).advance(14,'month');
  var date_range = ee.DateRange(start,end);
  var S1 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
    .filterDate(date_range)
    .sort('CLOUD_COVER')
    .filterBounds(p1).first()

  var S2 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
    .filterDate(date_range)
    .sort('CLOUD_COVER')
    .filterBounds(p2).first()

  var mosaic = ee.ImageCollection([ee.Image(S1), ee.Image(S2)]).mosaic();
  var image = mosaic.clip(TSL);

  return mosaic
};

var list_of_images = dates.map(fnc);
print('list_of_images', list_of_images);
var mt = ee.ImageCollection(list_of_images);
print(mt);
Map.addLayer(mt, {}, 'mt');

ERROR:
        List (Error)
        ImageCollection.fromImages: Attempt to create an ImageCollection with non-image elements.


texts:
What is the best way of geocoding IP addresses?

I have ~ 10,000 - 30,000 addresses to geocode once. The IP addresses are from all over the world and I only want the Long/Lat. Accuracy to the city level is not essential, but ideally it would be within 10-20 km.

Searching around I saw MaxMind which seems like a good option, but before I download and frequent myself with their API, are there any other recommended alternatives out there?


texts: Installing plugin from URL in XML doesn't move files in plugin folder

texts:
When I try to install a plugin using a repo I created the downloaded zip file doesn't place the files within the plugins folder correctly. They are placed outside the folder.
For example QGIS is expecting folder/files but when the .zip file is unzipped the files are being placed directory in the /plugins directory and not inside its own directory. I am thinking this is a QGIS 3.8 issue? Is there a different format of XML to use? Below is what I am currently using:

<?xml version = '1.0' encoding = 'UTF-8'?>
<?xml-stylesheet type="text/xsl" href="" ?>
<plugins>
    <pyqgis_plugin name="lll" version="1.0" plugin_id="227">
        <description>
            <![CDATA[QGIS tool - For common GIS tasks]]>
        </description>
        <about></about>
        <version>1.0</version>
        <qgis_minimum_version>3.0</qgis_minimum_version>
        <qgis_maximum_version>3.99</qgis_maximum_version>

        <file_name>networktool.zip</file_name>
        <icon>
        </icon>
        <author_name><![CDATA[prtdonn]]></author_name>
        <download_url>http://my-site.com/somePlugin.zip</download_url>
        <uploaded_by><![CDATA[jll]]></uploaded_by>
        <create_date>2015-05-16T15:31:19.824333</create_date>
        <update_date>2015-05-15T15:31:19.824333</update_date>
        <experimental>False</experimental>
        <deprecated>False</deprecated>
        <tracker>
        </tracker>
        <repository>
        </repository>
        <tags>
            <![CDATA[development,debugging,tools,LTE]]></tags>
        <downloads>20</downloads>
        <average_vote>5</average_vote>
        <rating_votes>10</rating_votes>
    </pyqgis_plugin>
</plugins>


texts: What software do you use as complementary to your GIS desktop/server?

texts: What software do you use as complementary to your GIS desktop/server in your work as GIS programmers/analysts and what programs are an absolute must you can't do without?


texts:
I am trying to use the tool to make the join a Mineral Database (excel) with my attribute table.  When I check the attribute table, it shows up as "NULL".  Why does it causing that and what did I do wrong ?

Can anyone please offer me any suggestions or ideas how do I fix it ? excel/database.


texts: Multiple-machine deployment scenario. Solutions?

texts:
I have a virtual machine (let's call it machine no.1) on which I have installed Microsoft SQL Express 2008 R+ ArcSDE 10.0+ ArcGIS Desktop 10.0 +ArcGIS Server 10.0. Everything is connected and works just fine.

Now I've installed ArcGIS Desktop 10.1 on another machine (machine no.2) and I want to connect to the Database Server on machine no.1. What exactly should I do? Install ArcSDE 10.1 on machine no.1?

If I try to Add a NEW Database Server on machine no.2, I get the error " Failed to connect to Database Server".

Thank you,


texts: I have an Normalized Difference Vegetation Index(NDVI) and Enhanced vegetation index (EVI) layer and I want to re-scale their data ranges to value ranges from 0–to–1, in order to optimally stretch the contrast (information contained in the image). I have analysed the histogram of the distribution of grey values and need to omit the lowest 0.5% and the highest 0.5% of the data, however I am not too sure how to go about doing this. I tried doing this in excel (manually deleting 0.5% min and 0.5% max values) and then completing the step in ERDAS model builder. The formula used in model builder (NDVI- new min)/(new max-new min), however my results were unsuccessful.


texts: military grid refrence system geocode reverse

texts:
I have lat/long and I get geocode like 34TCK9320739642 from mgrs geocode reverse in arcmap

how can convert that to address (country.city.state)


texts:
Doing some reverse geocoding for an application I am developping, but my qps seems capped at 1 qps.

Is this a hard limit, or is there a way to circumvent that?

I can't find any information on Mapzen.


texts: Convert Lambert Conformal Conic to WGS84

texts:
I am trying to convert my "Lambert Conformal Conic map with 2 standard parallels" to the basic/standard WGS84 projection that Google Maps use.

I have 30 of these maps which I need to convert. They are all in .tiff format with a corresponding .tfw file each.

What my final goal is:

Convert Maps to WGS84
Merge all 30 maps
Change them to tiles

For now I am only struggling with the conversion of projection.

Does anybody now of any ways/programs/techniques or anything (preferably free) to accomplish this?

I have heard/read about GDAL etc. but I do run Win Vista, not MAC.


texts:
Is there any command line tool to create a 3D terrain from a DEM SRTM HGT file like this?

I mean, given a camera position, a HGT file and a satelite image.

EDIT 1

First try to adapt mdsumner's answer:
This is the full code:

First, install all I need (Docker image from R-Base):

apt-get update & apt-get install -y libproj-dev xorg libx11-dev libglu1-mesa-dev libfreetype6-dev libgdal-dev --fix-missing

Install required packages:

install.packages(c("proj4"))
install.packages(c("reproj"))
install.packages(c("remotes", "quadmesh", "rgl"))
install.packages(c("raster"))
install.packages(c("rgdal"))

Load the SRTM DEM HGT file:

library(rgdal)
library(raster)
dem <- raster("/home/docker/S14W042.hgt")

Load a texture image and execute. My image is just a PNG file (must be a sub product of the HGT file - anything that I can extract using GDAL - or an ordinary image):

im <- brick("/home/docker/S14W042.map.png" )
qm <- quadmesh::quadmesh(dem, texture = im)

Error:

Error in quadmesh.BasicRaster(dem, texture = im) :
  texture must be a 3-layer raster with RGB values (in 0-255)

Please help here. NOTE: I can't auto-download anything. All files MUST be here when creating the image.

EDIT 2

I realized that I can use any X,Y,Z map server instead MapBox to get the imagery, including my own server.

First I try with ArcGIS online:

Give any API-KEY. Don't care.
Sys.setenv(MAPBOX_API_KEY="fdfsdfsdfsdfsdf")

u <- "https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/tile/{zoom}/{y}/{x}"

im <- cc_location(cbind(147, -42), base_url = u)

Done.


texts: Converting TIN to Multipolygon in PostGIS?

texts:
ST_Tesselate returns TIN geometry. But then I need to use other software to process geometry data and it only accept (multi-)linestrings and (multi-)polygons.

Is there way to convert TIN triangles to multipolygon?

TIN is 2-dimentional, no Z coordinate.


texts: Joinn shapefile attribute ID to point in R

texts:
I have a CSV file (about 1,000,000 georeferenced points) which contains point_id/latitude/longitude. I also have a shapefile which represents a grid cell_id/xmin/xmax/ymin/ymax.
I am trying to join the grid attributes (specifically the cell_id) to the points from my CSV file. Basically I want to know in which cell the points fall.

So my output would look like that point_id/latitude/longitude/cell_id.

Here is the code I have been working on based on this post and that post:

library(sp)
points<-read.csv(file = "/Users/input.csv", stringsAsFactors=FALSE)
coordinates(points)<- ~longitude+latitude
library(sf)
grid <- read_sf("/Users/grid.shp")
a.data <- raster::intersect(points,grid["cell_id",,,])
points$cell_id <- a.data$cell_id
write.csv(points, file = "/Users/points.csv")

However this does not gives me the cell_id but instead returns the point_id twice.

I have also tried using over but it yielded an error message.


texts: How can I import AutoCAD 3D solids into ArcGIS?

texts:
I have an AutoCAD drawing that contains contours and 3D solids.  I am able to import the contours (which i don't need) but not the 3D solids (which I do need).  I understand ArcGIS does not support AutoCAD 3D solids but is there a workaround workflow to do this?

EDIT:
It was pointed out that this type of file should be supported in version 9.3.  I can confirm that ESRI Help for 10.1 states that polygon and multipatch solid and 3D face is supported "on paper" in version 10.1 as well.

However in the dxf I have solids are ignored in ArcGIS (ArcScene or imported) but it can be opened in AutoCAD, DWG TrueView 2013, SketchUp 8.

DWG TrueView 2013:

ArcScene 10.1:


texts: Geofencing (Leave/Enter geofence Timing) Issue?

texts:
I already search SE already if the question has been raised before but I can't find one so I'll ask here anyway. I'm creating a small application which involve geofencing my GPS device doesn't have a geofencing capability on it's own which I hope it has this feature for accuracy and ease of programming.

I have no choice to do it on the server level, my question is since my device is sending NMEA/GPRMC on a two minute interval is there an algorithm that we can use to detect the exact timing the vehicle leave or enter the geofence?

my current code can detect the vehicle inside/outside the geofence area no matter what type of polygon it's just the timing issue I can't solve.


texts:
I am trying to construct a legend based on size using Carto VL, and am struggling to implement the solution posed in this Github thread.

https://github.com/CartoDB/carto-vl/issues/1260

I built a box to house the legend earlier in the code, here is the code I have so far that constructs the legend itself. What is this code missing?

const Legend = piceaLayer.viz.width.getLegendData(
                  width : ramp(($picea), [1, 2, 3]);
                 [{key: '10%', value: 10,
                   key: '50%', value: 50,
                   key: '100%', value: 100}]

               document.getElementById('legend').innerHTML = LegendList;


texts:
I'm working on some software that displays a map, like this:

![enter image description here][1]

I also have some export functions which exports the about image as .bmp/.Tiff/.Jpg/.PNG/.Wmf/.SVG

This is .TIFF export:

![enter image description here][2]

Now I also have to export it as GeoTiff. The software I am working on is in Delphi 7.

information

I have come across this steps to creating a geotiff, but since I am new to GeoTIFF, I am stuck after the first point.

question

Is it possible to convert the .TIFF file to .geotiff?

what I have tried

remotesensing geotiff

how to write geotiff files in delphi possibly-compressed

but I can't figure out how to create the meta data for the tiff to create the GeoTIFF


texts: Loading 'mod_spatialite' for SQLite in Jupyter Notebook (Windows)

texts:
I am trying to load mod_spatialite dynamic library for SQLite in Jupyter Notebook
import sqlite3

conn = sqlite3.connect('rl.db')
conn.enable_load_extension(True)
conn.execute('SELECT load_extension("mod_spatialite")')
conn.execute('SELECT InitSpatialMetaData(1);')

and I get the error
---------------------------------------------------------------------------
OperationalError                          Traceback (most recent call last)
<ipython-input-2-589a9312e5ee> in <module>
     14 conn.execute('SELECT load_extension("mod_spatialite")')
     15 conn.execute('SELECT InitSpatialMetaData(1);')
---> 16 conn.execute('SELECT load_extension("libspatialite")')
     17 conn.execute('SELECT InitSpatialMetaData();')
     18

 OperationalError: The specified module could not be found

I have tried the last solution from
https://stackoverflow.com/questions/8590599/sqlite-load-extension-fail-for-spatialite-in-python
but simply it doesn't work, I got even more complicated errors related to other libraries dependent on sqlite3.dll.
Does anyone succeed in working with spatialite extension via Jupyter Notebook and knows how to properly load mod_spatialite?
Python v3.9.1 (tags/v3.9.1:1e5d33e, Dec  7 2020, 17:08:21) [MSC v.1927 64 bit (AMD64)]

SQLite v3.33.0


texts: I am new user to Geomajas. I have created Geomajas plugin and I am using geotools to add shapefiles. and I want to create buffer on layer(road layer, or polygon layer), but I am not getting exact steps to create buffer on Geomajas map, so please help me to to  add buffer to map.


texts: Accessing watershed boundary for lakes in R

texts:
I have had success accessing unique database IDs and polygons for various lakes in the US. I'm now looking to get polygons for the watersheds associated with these lakes, but have stumbled a bit.

Example- I get a lake polygon from a GPS point using the Hydrolinks package:

library(hydrolinks)

Lat<-37.70127
Lon<- -97.0508

Lake <-link_to_waterbodies(Lat, Lon, "Name", dataset = c("nhdh", "hydrolakes", "nhdplusv2"), buffer = 50)

The result is Santa Fe Lake in Kansas. (ID = 45827774)

Using the "nhdR" package, I can get the polygon associated with this lake

library(nhdR)
dt <- nhd_load(c("KS"), c("NHDWaterbody"))
Poly<-dt[dt$Permanent_Identifier == 45827774]

I'm using results from the streamstats API as my "true" watershed (see screen shot from streamstats API for this lake).

The reachcode for the "Poly" object is "11030017001699". And if I try to calculating watershed data from the first 10 digits of the reachcode, the basin doesn't match what I get from the streamstatas API.

I tried using the pour point for this lake (Hydrolakes lists this as -97.05319   37.709), but calculating a watershed based on these points only gets me a small section of the true watershed.

library(streamstats)

#Pour Points
PLon<- -97.05319
PLat<- 37.709
ws1 <- delineateWatershed(xlocation = PLon, ylocation = PLat, crs = 4326,
                      includeparameters = "true")

So I'm assuming that the accuracy of these pour points is not good enough to use streamstats to calculate entire lake watersheds.

Using tools from nhdplusTools package works for this lake

library(leaflet)
library(nhdplusToos)

 #Step 1
  point <- sf::st_sfc(sf::st_point(c(PLon, PLat)), crs = 4326)
  #Get nhdplus_id from point
  ID<-discover_nhdplus_id(point)

  #Get Basin
  nldi <- list(featureSource = "comid", featureID = ID)
  basin <- get_nldi_basin(nldi_feature = nldi)

  leaflet() %>%
    addPolygons(data = basin) %>%
    addCircleMarkers(data = point,
                     #lat = ~lat, lng = ~long,
                     color = "red",
                     radius = 2) %>%
    addTiles() #Add tiles adds defualt map

Which seems to work well. But I don' always have pour points for lakes, or some of my pour points do not give me a full watershed. Using random GPS points from within a lake polygon gives me unreliable results.

 Lat1<- 37.72797
 Lon1<- -97.82922
 ID = 135039628

Plugging the above coordinates into the "get_nldi_basin" function only gets a small aspect of the entire watershed

And if I used these coords for a small pond:

Lat2<- 39.27623
Lon2<- -94.89211
ID = 159110681

I get a huge watershed

So using the "get_nldi_basin" seems to be limited to my ability to get an accurate pour point. I'd really prefer a method that gives me an accurate watershed boundary for an entire lake using the entire lake polygon, or a match to the permanent NHD ID.


texts: How to gdal_translate a NOAA GFS grib2 file and cut a piece out with -projwin (coordinate system chaos)?

texts:
So I'm trying to convert NOAA GFS gribs files to PNGs with gdal_translate. My grib2s have a pretty big area undefined with wgrib, which appears in white on the resulting image. I have been trying for hours to get the crop -projwin right but the coordinate system chaos is melting my brain. Those grib2 files come with their own coordinate system which gdalinfo shows as
    Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: gfs.t00z.pgrb2
Size is 1440, 721
Coordinate System is:
GEOGCRS["Coordinate System imported from GRIB file",
    DATUM["unnamed",
        ELLIPSOID["Sphere",6371229,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (-0.125000000000000,90.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Corner Coordinates:
Upper Left  (  -0.1250000,  90.1250000) (  0d 7'30.00"W, 90d 7'30.00"N)
Lower Left  (  -0.1250000, -90.1250000) (  0d 7'30.00"W, 90d 7'30.00"S)
Upper Right (     359.875,      90.125) (359d52'30.00"E, 90d 7'30.00"N)
Lower Right (     359.875,     -90.125) (359d52'30.00"E, 90d 7'30.00"S)
Center      ( 179.8750000,   0.0000000) (179d52'30.00"E,  0d 0' 0.01"N)

I basically want to crop to an area from lon -40 to 80 and lat 30 to 90 and then output a file in EPSG:3857.
My script looks something like this (minus the -projwin):
${GDALXLATE} -of png -b 128 /grib/gfs.t06z.pgrb2 /out/test.png
${GDALWARP} -overwrite -setci -t_srs EPSG:3857 /out/test.png /out/test.tif
${GDALXLATE} -of png /out/test.tif /out/test.png

I tried a whole bunch of different -projwin parameters but none of them produced results even close to what I'm looking for... Maybe somebody who has little coordinates running through their bloodstream could help me out here? I'm kinda desperate...
-Stefan


texts: How do you export to georeferenced PDF from QGIS 2.0?

texts:
When I use the print composer in QGIS, there is no option to export to georeferenced PDF using GDAL.  Do I need to export from somewhere else or how do I go about exporting to georeferenced PDF?  I am using the newest version of QGIS.

Related questions for older versions: Is it possible to export geospatial PDF from QGIS via GDAL/OGR?


texts: Values in symbology tab of layer properties not showing up

texts:
I am trying to colorize (quantify) a chloropleth map based on rent price data in ArcMap 10.2.1. The rent price value is from my .xlsx excel table, which is joined with my shapefile. However, the rent price value is missing from the drop-down list of values under the symbology tab. I know that the data was joined without any errors, so I'm not sure why I'm not able to color-code my map based on the rent price data.

Does anyone have an idea what I am missing?

EDIT: The rent price was listed as a string value and so I added a new numeric field to the attribute table of the census tract shapefile (shown below), but I'm not sure how to get the values from the string column into the numeric column.


texts: Constraining point geom OL3 feature marker movement to that of a multi line string

texts:
I have a vector layer in OL3 that I generate using a .gpx file.  I have created another vector layer which includes 2 point geometry features.

I want to be able to click and move the marker features, but constrain their movement so that they can only be moved along the multi point line string of the associated gpx feature vector layer.

I was hoping that the Snap interaction would work, but this messes up the dragging of the map itself.

The code below shows how I am doing this.  The 2 marker features are created and added to the same OL source that contains ALL of the features from the GPX file I loaded.

        //create the start and end marker features
    var startMarkerFeature  = _createMarkerFeature(startIconCoords, startIconStyle, "Start Marker");
    var endMarkerFeature    = _createMarkerFeature(finishIconCoords, finishIconStyle, "End Marker");

    var markerFeatures = new ol.Collection([startMarkerFeature, endMarkerFeature]);

    //Add the start and end marker features to the track point layer containing all of the GPX track features
    gpxTrkPointLayer.getSource().addFeatures(markerFeatures);

    //Make the 2 markers draggable
    var markerDragInteraction = new ol.interaction.Modify({features: markerFeatures, style: null});

    //Callback for when user has stopped dragging the marker
    markerDragInteraction.on('modifyend', markerDragged);

    map.addInteraction(markerDragInteraction);
    map.addInteraction(new ol.interaction.Snap({source: gpxTrkPointLayer.getSource()}));

This code enables me to drag the 2 markers around, but the drag movement of the 2 features is not locked to the features contained in the same source.  I know I am missing something here, but cannot deduce what.


texts: Land use change evaluation using Molusce in QGIS?

texts:
When i upload my two raster files in Molusce (initial and final one) i get the message that the files have different geometry.
Their coordiante system is the same as well as their tif size.
In order to get these raster files, I have converted two vector files.
Vector file for the initial year is made out of classified polygons which have been done using an old aerial photo. Vector file for the final year is made out of classified polygons which have been done using google earth sattelite image.

Is it possible that the problem is slightly different shapes of the polygonized area?


texts:
What I need is a function for Non Earth Coordinates in QGIS (both projects and layers). They should be mappable to any projected coordinate system (like SWEREF99TM).

Why: In the Mining industry we often use Non Earth systems for different reasons. Like Geophysics we want to define North in a convenient direction. Could be South East we don't care about the real world at this stage. The same if an ore body is going south-east we may define that direction as north.

In many softwares (Microstation or MapInfo) we use we can have on the fly conversion between lets say myCoordinateSystem1 (cartesian system) to EPSG3006 no problem. But I have not found a good function for this in QGIS.

I wonder if the Proj4 limits QGIS so you can't do this?
Any solutions for this?