texts:
I am new to gdal and have a simple question. I want to transform an image from OpenStreetMap projection to Nasa Blue Marble projection. What I understand is that OpenStreetMap projection is EPSG:3857 and Blue Marble projection is EPSG:4326 (http://www.paulillsley.com/gia/index.html).

So, I simply downloaded a sample 256*256 tile like this: http://a.tah.openstreetmap.org/Tiles/tile/0/0/0.png and applied the following command:

gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 "0.png" "o_bluemarble.png"

I received the following error:
"Copying color table from 0.png to new file.
ERROR 1: Unable to compute a transformation between pixel/line
and georeferenced coordinates for 0.png.
There is no affine transformation and no GCPs."

From where should I start?

Update 1

I got it. I need to define the bounds. Firstly, I got the bounds:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
-180 85.05
-20037508.3427892 20036051.9193368 0

I used 85.05 as it is the approximation to infinity in Spherical Mercator used by Bing. Then, I used the following commands:

./gdal_translate -of Gtiff -co "tfw=yes"  -a_ullr -20037508.3427892 20036051.9193368 20037508.3427892 -20036051.9193368 -a_srs "EPSG:3857"  "input.png" "input_tfw.tiff"

./gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 -ts 256 128 "input_tfw.tiff" "output.tif"

The problem now is that by comparing an image downloaded from openstreetmap with Blue Marble, the poles are not aligned properly. I mean here Greenland, not the far pole. What do think the problem is?


texts: Adding filter functions in place of a literal expression in GeoTools?

texts:
I'm writing some code to use filter functions in geotools.

(attName >= replacementFunc() AND title LIKE replacementFunc()) AND attName2 = replacementFunc()

In the area that says LIKE I get a parsing error. If I change it to EQ there are no problems. If I wrap it in quotes, then I get a literal 'replacementFunc()' inside the output. What are my options?


texts: Select a specific line length after and before a point

texts:
As you can see from the attached image, I have a shapefile of river network (line) and a point shapefile (this one represents the middle point of transects carried out in the area). I'm trying to select 300 m of the lines, both before and after the points. The aim is to finally have a 600 m linear buffer in which each point is in the center.

I've already tried with a 300 m buffer for each points, but that doesn't work because some river segments (inside the buffer) are longer/shorter than 600 m.


texts:
I have a custom python tool which creates a number of rasters and writes them to C:/Documents/ScriptTool folder as tif files, and also adds them as layers to the active map.

At the startup of the tool, I remove all layers from the map via:

aprx = arcpy.mp.ArcGISProject("CURRENT")
activeMap = aprx.activeMap
layers = activeMap.listLayers()

for layer in layers:
    if not "Map" in layer.longName:
        activeMap.removeLayer(layer)
aprx.save()

And I attempt to delete all files which are in the C:/Documents/ScriptTool folder, because I do not want them to linger for memory concerns and also because ArcGIS likes to crash if you try to do so.

try:
    shutil.rmtree("C:/Documents")
except:
    arcpy.AddMessage("Failed to delete all files.")

However, the .tif files themselves are not being deleted (Though the .tif.aux files and the .tfw files are). When I attempt to manually delete the files through the file explorer, I'm given the message: The action can't be completed because the file is open in ArcGIS Pro., despite having removed those layers from my active map.

If I close and re-open ArcGIS Pro, I'm able to delete the files and run my tool once again, but obviously this is non-ideal. Can I release the hold ArcGIS Pro has on these files so they may be deleted, programmatically?

I have also tried using arcpy.Delete_management() to remove the layers, with identical results.


texts: Creating cloud free imagery and applying NDVI calculation to NAIP, Sentinel 2, and L8 collections using Google Earth Engine?

texts:
I'm trying to get cloud and shadow-free (or as close to cloud free) scenes for Landsat 8, NAIP, and Sentinel 2 (specifically sentinel 2 10 meter resolution), and then apply the NDVI calculation to each set.

After that, I'd like to look at those values within specific polygon geometries I've created, but that's a separate thing I'd like to try on my own first!

I'm fairly new to GEE, and though I've watched the tutorials I still have some trouble with basic things it seems.

Here is just the portion I have so far for Landsat 8, though it's not running error-free.

// add three imagery layers to neighborhoods and filter date
var l8filtered = landsat8.filterDate('2016-10-01', '2016-10-31');
                         .sort('CLOUD_COVER');
                         .first());

// send to console/add to map
print(l8filtered, "l8 filter")
Map.addLayer(l8filtered,{},"L8")

//create NDVI image
var l8_ndvi = l8filtered.normalizedDifference(['B5','B4']);
Map.addLayer(l8_ndvi,{min:-0.1, max:0.1, palette:['FFFFFF', '339900']},'L8 NDVI');

For the cloud shadow - I've been sleuthing a did come across this post, however it's not quite making sense to me.

For the Landsat 8 imagery specifically - am I right in defining the bands in the l8_ndvi variable, or should I be calling the bands out before that?


texts: Can I use RGeo::Point data in PostGIS queries without multiple conversion?

texts:
How can I represent RGeo::Point ( CAPI implementation ) in some way where I don't have to convert it once more in Sql?

Now I have something like this:

class Poly
  scope :enclosing ->(geom) do
    where ["ST_Contains(poly.geom, ST_GeomFromText( ? , 4326))",
           geom.as_text]
  end
end
...
p = RGeo::Geos.factory(srid: 4326).point(10, 20)
Poly.enclosing(p)

This way the point first gets represented as "POINT (10.0 20.0)" , then converted to geometry by ST_GeomFromText. Is there a shortcut? PostGIS geometry items look like 0101000020E6100000CC79C6BE64E13D4046F6F708EB8F4E40 and have SRID already encoded, but RGeo::Point.as_binary returns something very different and not acceptable by PostGIS.


texts: OpenLayers - edit properties of GeoJSON vector source

texts:
I have a GeoJSON vector file with some polygons representing country borders. I have added a property called "processed" to each country which defaults to 0 (i.e. false). I want to create a function that changes this value to 1 (true) when i select a country on OpenLayers map and click a button. And i want to overwrite old data with this edited data or make new data exportable to server-side.

Vector data source looks like this:
https://raw.githubusercontent.com/mleontenko/openlayers-edit-properties/master/borders1.geojson

Full html file is here:
https://github.com/mleontenko/openlayers-edit-properties/blob/master/map1.html

Script code (commented missing parts) looks like this:

  var vectorSource = new ol.source.Vector({
    url: 'https://raw.githubusercontent.com/mleontenko/openlayers-edit-properties/master/borders1.geojson',
    format: new ol.format.GeoJSON()
  });

  var map = new ol.Map({
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      }),
      new ol.layer.Vector({
        source: vectorSource
      })
    ],
    target: 'map',
    view: new ol.View({
      center: ol.proj.fromLonLat([4.8, 47.75]),
      zoom: 5
    })
  });

  // a normal select interaction to handle click
  var select = new ol.interaction.Select();
  map.addInteraction(select);

//button that changes 0 to 1
var PushButton = document.getElementById('store-data');
PushButton.addEventListener('click', function(){
  if(select){
    // missing code that changes "processed" value from 0 to 1

  }else {
     window.alert("You have not selected anything");
    }
 });

So, how can i change value of "processed" property from 0 to 1 and store changes to GeoJSON to "update" the map?


texts: Creating Data Clock in ArcGIS Tracking Analyst?

texts:
I am attempting to create a data clock, under the "Tracking Analyst" toolbar.  I have created a tracking layer, but for some reason the "Create Data Clock" option is grayed out.

Does anyone know why this is?


texts: Adding ArcGIS Online Services to QGIS using MAC OSX?

texts:
I just read Adding ArcGIS Online Services to QGIS?, but I do not have the options mentioned.

These are my options under Add Layer... Can't seem to find Add ArcGIS Map Server Layer. I'm running 2.18 on a Mac. There must be a difference between mac and pc versions? Is there a way to add ArcGIS Map Server Layer on a Mac??


texts: Reprojecting in GeoServer with custom transformation parameters

texts:
I already set up a raster layer in GeoServer, store type is GeoTIFF with the EPSG:2462 and I'm trying to use it as WMS in AutoCAD and QGIS using the target EPSG:32634.

If I insert the same raster as a file based it is reprojected as it should be, but if I insert it as WMS through GeoServer I get it shifted several meters east.

I see from the QGIS coordinate system that it already has the transformation parameters correctly for that EPSG but I am not sure that when you insert a WMS raster layer it is reprojected by QGIS or through the GeoServer?

If it is reprojected by GeoServer how can I manage it or how can I use only the QGIS reprojection instead so everything to be OK as it is when you use a file based raster?

This is the CRS in QGIS:
+proj=tmerc +lat_0=0 +lon_0=21 +k=1 +x_0=4500000 +y_0=0 +ellps=krass +towgs84=-44.183,-0.58,-38.489,2.3867,2.7072,-3.5196,-8.2703 +units=m +no_defs


texts: FID attribute for each circular buffer polygon

texts:
I have a point shapefile "Negative_vel". With the command Analysis tools-Proximity-Buffer and selecting 50m as a linear unit and dissolve type ALL I created a new .shp "Negative_vel_buffer50m". But this new shapefile (with the buffer circular polygons) has no attribute table.
I need an FID for each circular buffer polygon so that to check with Join Points to Layer how many points from the initial shapefile "Negative_vel" are within each buffer polygon. Because I need to keep only the polygon buffers with more than 3 points inside them.
It is not of my interest to keep buffer polygons with just one point inside.

So, I need a buffer .shp dissolved with attribute table so later on to Join and count the points inside each polygon.


texts: How to convert the file with projections in meters to Lat-Long in QGIS ar Erdas Imagine?

texts:
I have file in Lambert Conformal conic format for the data of Australia, Its coordinates are in meters. I want to convert it to Latitude-Longitude format. Which Projection, Spheroid and Datum I should be using if I want to convert the projections in QGIS or Erdas Imagine?

The current Specification of the file is as follow:
Projection: Lambert Conformal Conic
Spheroid: GRS 1980
Datum: GDA94-ICSM
Map Units: Meters

Many Thanks in Advance for the Answer.


texts: I have point with an attribute for each point's 0-360 direction. I need to render these points in Leaflet with arrows pointing in their respective directions. Is there a plugin or code I can use to accomplish this?


texts: How to remove geoprocessing history in ArcGIS 10.2.2 for Desktop?

texts:
I have grown dependent upon the awesome add-in called Clear Geoprocessing History.

It does not appear to work with my upgrade to 10.2.2

Does someone know of a different route?

Or perhaps how to modify the tool to make it work?


texts: Extracting a single grid from a mdb file?

texts:
I'm trying to open ISRIC-WISE derived soil properties (available as MDB file from http://www.isric.org/data/) as a layer in ArcMap without much success. I know MDB are MS Access files that can be read in ArcMap as a geodatabase, but I'm not familiar with MDB files in GIS. I got to open the MDB file and see each one of the 20 derived soil properties as a spreadsheet, but not as a layer.

1) How can I see each of the 20 derived soil properties as a layer in ArcMap 10?

2) Once imported the 20 layers, how can I extract values of some X & Y points from each one of the 20 layers at the same time?


texts: Force to Display All the Trail Names on a Map in QGIS

texts:
I have been testing ways to force all the MTB trail names to be visible on an OSM map I am styling in QGIS.

I can see them as red unplaced labels and there is room on the map to display them yet they do not appear.

I have set in the Placement Engine to show all labels. Played with the overrun feature, priority and discourage from covering settings yet nothing works fully. The font is small enough to certainly allow other names to show with ample space around them.

Does anyone know how to make ALL these hidden red trail names appear?

Dan


texts: ArcObjects Memory Leak with IWorkspaceEdit?

texts:
I have ArcGIS 10.1 installed.
I have created a function for copying all of the features from one shapefile to another using IFeatureBuffer, now every time the function when it is called, it takes some memory, but doesn't release that much, meaning on every function call a memory leak happens. What i have observed is that the memory leak is dependent upon the size of the shapefile, I imported the shapefile with around 34000 features, and the memory leak was much greater than that with a smaller shapefile. Following is the minimal code version of this behaviour:

public void CopyFeatures(string srcFile, string srcDir, string outFile, string outDir)
{
    IWorkspaceFactory workspaceFactory = new ESRI.ArcGIS.DataSourcesFile.ShapefileWorkspaceFactory();
    IWorkspace workspace = workspaceFactory.OpenFromFile(srcDir, 0);
    IFeatureClass featureClass = ((IFeatureWorkspace)workspace).OpenFeatureClass(srcFile);
    IFeatureCursor featureCursor = featureClass.Search(null, false);
    IFeature feature = featureCursor.NextFeature();

    IWorkspace outputWorkspace = workspaceFactory.OpenFromFile(outDir, 0);
    IFeatureClass outputFeatureClass = ((IFeatureWorkspace)outputWorkspace).CreateFeatureClass(outFile, featureClass.Fields, null, null, esriFeatureType.esriFTSimple, "Shape", "");

    IWorkspaceEdit workspaceEdit = ((IDataset)outputFeatureClass).Workspace as IWorkspaceEdit;
    if (!workspaceEdit.IsBeingEdited())
    {
        workspaceEdit.StartEditing(false);
        workspaceEdit.StartEditOperation();
    }

    IFeatureBuffer outputFeatureBuffer = outputFeatureClass.CreateFeatureBuffer();
    IFeatureCursor outputFeatureCursor = outputFeatureClass.Insert(true);

    while (feature != null)
    {
        for (int i = 0; i < outputFeatureBuffer.Fields.FieldCount; i++)
        {
            IField field = feature.Fields.Field[i];
            string myName = field.Name;
            if (field.Editable == true)
                outputFeatureBuffer.Value[outputFeatureBuffer.Fields.FindField(myName)] = feature.Value[feature.Fields.FindField(myName)];
        }

        outputFeatureCursor.InsertFeature(outputFeatureBuffer);
        ReleaseComReferences(feature);
        feature = featureCursor.NextFeature();
    }
    outputFeatureCursor.Flush();
    if (!ReferenceEquals(workspaceEdit, null))
    {
        workspaceEdit.StopEditOperation();
        workspaceEdit.StopEditing(true); // This line causing the increase in memory but the same memory doesn't get released.
        ReleaseComReferences(workspaceEdit);
    }

    ReleaseComReferences(featureCursor);
    ReleaseComReferences(workspace);
    ReleaseComReferences(featureClass);
    ReleaseComReferences(outputFeatureCursor);
    ReleaseComReferences(outputFeatureBuffer);
    ReleaseComReferences(outputWorkspace);
    ReleaseComReferences(outputFeatureClass);
    ReleaseComReferences(workspaceFactory);
}

public void ReleaseComReferences(object refernce)
{
    if (refernce != null)
    {
        int refsLeft = 0;
        do
        {
           refsLeft = System.Runtime.InteropServices.Marshal.ReleaseComObject(refernce);
        }
        while (refsLeft > 0);
    }
}

Now, each time the CopyFeatures function is called, some memory doesn't get released, and the memory that is being used keeps on increasing. So any suggestion, how i can fix this memory leak issue? I have searched a lot on the internet but wasn't able to find a viable solution, every solution that i tried wasn't working, and still some memory wasn't being released.

Also i have checked with the copyfeatures tool and the memory issue doesn't occur, but the above code is just an excerpt from my original code, I also need to perform some logic and some checks on every feature, due to which I cannot use copyfeatures in this scenario.


texts:
I have 2 Sentinel files that I join in R:

library(raster)
y = raster("./L1C_T36MZE_A016248_20180802T075248.tif")
z = raster("./L1C_T36MZE_A016291_20180805T080721.tif")

ss = raster::merge(y,z)
plot(ss)

This gives a file with a seam (the diagonal line) where the images have spurious edge values.

What is the best way to deal these values?
I have tried:

sp = click(a, xy = T)
spp = spPolygons(as.matrix(sp)[, 1:2])
plot(crop(a, spp))

but this did not help.


texts: QGIS attributetable: show $x and $y in other CRS than layer

texts:
I have a point file with an attribute table. In the attribute table I am able to display the X and Y coordinates using the field calculator. These coordinates are in ETRS89. Because the end-user of the data uses a different system (EPSG:28992) I would like to display two more columns in my attribute table with the X and Y in EPSG:28992.

I thought I'd use the transform function but I can't get it to function. I would like the values to change if I move the location of the point in my chart.

How do I fix this?


texts:
Is there a way to compute isochrones for public transport while using the plug-ins available in QGIS? Or is there another way for making them?

Note: I am not a software engineer which means that some of your answers might exceed my field of study regarding programing or handling specific software.


texts: How to find closest point without crossing obstacles?

texts: I have a set of points and I'd like to be able to find the nearest point from another layer. The problem is that I need to find the closest point without crossing an obstacle (such as a river, highway, etc). Any idea how to do this?


texts:
I very new to all the python-QGIS issue. I try to run this one line code in order to know if the PYTHON 2.7 recognize the modules of QGIS:

import qgis.core

and got this error:

Traceback (most recent call last):
  File "G:\desktop\test1.py", line 1, in <module>
import qgis.core
ImportError: No module named qgis.core
>>>

i red

http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html#python-plugins

Where is the python library of QGIS 2.8

but didn't understand what should i do?


texts: What coordinates: 6 digits, 8 decimals and how to convert to Lat & Lng

texts:
I have found a public GeoJSON dataset that uses coordinates I can't identify.
The first tuple is:

[146476.39635656,165501.27742543]

This tuple is located in Brussels, Belgium (UTM 31) and is close to [4.35659887, 50.83598669].

What format is it?
What is the latitude & longitude of the example above?


texts:
I am extremely new to PostGIS, and I am interested in visualizing a large GPS dataset. The dataset is updated daily and contains millions of points in the format (EPSG:4326):

Device ID    Timestamp    Latitude    Longitude

Mapbox Vector Tiles seem like a good approach for visualizing GPS tracks on a map, but I do not have good intuition as to the most efficient design. Currently, I have the following approach:

Use St_MakeLine to convert the GPS coordinates for each device to a line geometry, ordering by the timestamp
Simplify the route using ST_Simplify to reduce the number of points in each line
Transform each line into a Mapbox Vector tile using ST_AsMVT and ST_AsMVTGeom

This seems to work reasonably well with a small dataset, but does not scale well as the number of points for each device grows (the St_MakeLine and ST_Simplify calls become more and more expensive).

How is this typically handled? Should the dataset be further broken down into tracks, resulting in smaller line geometries? Should it be stored as lines instead of points? According to this question, adding new data is slow, but that may be acceptable with a large boost to the visualization performance.


texts: Identify countries within a bounding box using R

texts:
Using R and relevant GIS packages, I would like to determine which countries are within an arbitrary bounding box and extract the corresponding ISO codes. For example...

require(maps)
require(rgdal)
require(raster)

e1<-extent(c(-109,-79,10,28))
bbox1<-bbox(e1)

How can I then use this bbox or four coordinates to determine which countries recognized in the GADM database fall within this bounding box?


texts: Replace OpenLayers draw form select options by button

texts:
I am trying to replace the openlayers draw features form by button. I have created three buttons for point, line and polygon features draw. When clicking the first button (value = point) it works fine but clicking second button it doesn't change the drawing value properties (like value = line).
Below is my code:

    <input id="type" type="button" value="Point">
    <input id="type" type="button" value="Linestring">
    <input id="type" type="button" value="Polygon">

    ...

    var typeSelect = document.getElementById('type');

    document.getElementById('type').onclick = function addInteraction() {
      var draw = new ol.interaction.Draw({
        source: source,
        type: typeSelect.value
      });
      map.addInteraction(draw);
    }

    /**
     * Handle change event.
     */
    typeSelect.onchange = function() {
      map.removeInteraction(draw);
      addInteraction();
    };

    addInteraction();


texts: How to save the raster file with the elevation data and convert to a asc/dem in arcgis pro?

texts:
Beginner arcgis pro user here...
I have subset and digitized a satellite imagery and I have the ground elevation data (first image). How do I save the data (integrated with the elevation) as a an asci/dem/grd
I need this as a surface file that I can use for Visual modflowflex(Second image)


texts: Model Builder/Python: batch raster calculation

texts:
I wish to carry out a raster calculation on a batch of rasters.
I want to add rasters by month, as identified by their file name.

i.e. all rasters within the month of January should add together and output a single jan.png/raster file. Then all the rasters within Feb should add together and produce a single FEB.png.

Usually I would use raster calculator, I would prefer to automate it via either python or modelbuilder.

My plan outline is as follows:

Input Files. 36 Weekly raster files
5jan15-11jan15.png
12jan15-18jan15.png
19jan15-25jan15.png
26jan15-31jan15.png
01feb15-07feb15.png....

Raster operation.
(raster1+raster2+raster3+raster4)/(raster1+raster2+raster3+raster4) = month.png

output. 12 rasters for each month.
jan.png
feb.png
march.png

is it possible to write a python expression or create a modelbuilder which will add rasters by month using their filename?

Current Attempts

I have built a model, however this only runs on the rasters I provide it, requiring me to write a new expression each time. In other words it's the same as manually running raster calculator and writing a new expression for each month, I wish to automate this via modelbuilder.


texts: If Statements in ArcPy Select Feature SQL Query

texts:
I'm trying to select pipes with specific codes. These codes are 8 digits long and I only need to select certain combinations. For example, pipes that have codes with the number 8 in the 7th position of the code need to be selected, meaning that a pipe with a code such as 02004580 would be selected. In total there are 4 different cases of codes that would require selection.

I don't know much about SQL but I tried to use the SelectLayerByAttribute tool and then input a query that would select pipes with codes that had an 8 in the 7th position but it did not work:

SELECT*FROM w_pipe WHERE:
SELECT SUBSTRING(CRITICAL_CODE,7,1) =8

The "SELECT*FROM w_pipe WHERE:" is what comes up automatically in the Select by Attributes window. I have a feeling the SelectLayerByAttribute tool might not support the SQL statement that I entered but I'm thinking I just didn't input the proper code.


texts: Using GDAL in Python class results in access violation (0xC0000005)

texts:
I am getting Access Violation error when running (some) GDAL functions from a class method, while running a straight script works fine.
Process finished with exit code -1073741819 (0xC0000005) at point of failure, which corresponds to Access Violation.
Code that fails:
from osgeo import gdal

class GeotiffReader:
    def __init__(self, file_path):
        self.file_path = file_path
        self.bands = []
        self.read_geotiff()

    def read_geotiff(self):
        ds = gdal.Open(self.file_path)  # works
        band_count = ds.RasterCount
        for i in range(1, band_count+1):
            self.bands.append(ds.GetRasterBand(i))  # this works too

    def construct_nn_input_data(self):
        arras = []
        print(sys.path)
        for band in self.bands:
            arr = band.ReadAsArray(0, 0, 500, 500)  # POINT OF FAILURE
            arras.append(arr)

        return arras

gtread = GeotiffReader("mapa_r_g_b_ir_mask.tif")  # or other multi-band tiff
gtread.construct_nn_input_data()

At the same time, when I simply define a method all works fine:
from osgeo import gdal

def read_arrs(file_path):
    bands = []
    ds = gdal.Open(file_path)
    band_count = ds.RasterCount
    for i in range(1, band_count+1):
        bands.append(ds.GetRasterBand(i))

    arras = []
    for band in bands:
        arr = band.ReadAsArray()
        arras.append(arr)

    return arras

read_arrs("mapa_r_g_b_ir_mask.tif")

I am sitting on a Windows 10 machine 64bit using Python 3.7.9 and GDAL 3.2.0 and am working on IntelliJ IDEA 2020.2.
Digging deep with debugger at assembler level showed the program really tries to execute some random part of memory.
I have tried working with PATH system variables, resorted removing both Python and GDAL, clearing registries, re-installing basically from scratch.
Most odd to me is that it works (with exact same settings, sys-variables, python, etc.) if we call ReadAsArray from a method, but I cannot wrap the code in a class.
What am I missing?


texts: Showing transparency slider by default in Layers Panel of QGIS?

texts: Is it possible to have a transparency slider shown by default in the Layer window, i.e. directly after adding a new layer, instead of having it to enable for each layer manually via layer-properties?


texts: How to show popup info after click to marker?

texts:
this is my code with some marker. I would like to show square with marker informations after click to some marker on map like there. I get piece of js code from this page, write to my code and tried it. But doesn´t work with info squares. That can you help me to choose correct part of code from page for show info square about marker - for example with text "Hello world!" after click to some marker? Thank´s

index.html

<!DOCTYPE html>

 <title>Working with Openlayers</title>
 <link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">

index.js (Just isn´t correctly running)

var baseMapLayer = new ol.layer.Tile({
  source: new ol.source.OSM()
});
var map = new ol.Map({
  target: 'map',
  layers: [ baseMapLayer],
  view: new ol.View({
      center: ol.proj.fromLonLat([-74.0061,40.712]),
          zoom: 12 //Initial Zoom Level
        })
});

var marker = new ol.Feature({
  geometry: new ol.geom.Point(
    ol.proj.fromLonLat([-74.006,40.7127])
     ),  // Cordinates of New York's City Hall
});

var marker2 = new ol.Feature({
  geometry: new ol.geom.Point(
    ol.proj.fromLonLat([40.7127, -74.006])
     ),  // Cordinates of New York's City Hall
});

marker.setStyle(new ol.style.Style({
  image: new ol.style.Icon( ({
    color: '#ffcd46',
    crossOrigin: 'anonymous',
    src: 'dot.png'
  }))
}));

marker2.setStyle(new ol.style.Style({
  image: new ol.style.Icon( ({
    color: '#cdff46',
    crossOrigin: 'anonymous',
    src: 'dot.png'
  }))
}));

//function for show marker info
      map.on('singleclick', function(evt) {
        var name = map.forEachFeatureAtPixel(evt.pixel, function(feature) {
          return 'Hello world';

        })

          container.style.display="block";
          var coordinate = evt.coordinate;
          content.innerHTML = "Hello";
          overlay.setPosition(coordinate);

      });
      map.on('pointermove', function(evt) {
        map.getTargetElement().style.cursor = map.hasFeatureAtPixel(evt.pixel) ? 'pointer' : '';
      });
//to there

var vectorSource = new ol.source.Vector({
  features: [marker, marker2]
});

var markerVectorLayer = new ol.layer.Vector({
  source: vectorSource,
});
map.addLayer(markerVectorLayer);


texts: Add ID if one feature is selected by another using PyQGIS

texts:
My inputs: line shapefile with features inside and polygon shapefile with other features inside.  What I am trying to achieve is to select from line feature every feature which is inside each polygon: if yes, it is added to the new column unique number.

So far I managed to create this script:

from PyQt4.QtCore import QVariant

cLayer = iface.mapCanvas().currentLayer()
bLayer = QgsMapLayerRegistry.instance().mapLayersByName("Mean coordinates")[0]

cLayer.dataProvider().addAttributes( [ QgsField("ID", QVariant.Int) ] )
cLayer.updateFields()
idx = cLayer.dataProvider().fieldNameIndex( "ID")
selection = []
n=1
for pt in cLayer.getFeatures():

    point = pt.geometry()
    for pol in bLayer.getFeatures():
        pol = pol.geometry()

        if pol.contains(point):
            cLayer.setAttribute(n, "ID")

        n+=1

It creates ID field, but my ID column is not updating. I am also not sure if my "select by location" query for each individual polygons works. Can you help me?


texts: Creating duplicate line set distance away from a polyline in QGIS?

texts:
I'm new to QGIS and not really sure what I'm doing; I have a polyline of the UK coastline and I want to add another line 30 km offshore.

Is this possible?


texts:
I have a script that iterates through layer subgroups in QGIS 3.6, opens a layout with the with the same name as the subgroup and exports a PDF. What I am struggling with is to close the active layout window after the export.
Any ideas or other ways of doing such task?

from qgis.core import QgsLayerTreeGroup, QgsLayerTreeLayer
import os,sys

root = QgsProject.instance().layerTreeRoot()
project = QgsProject.instance()
projectLayoutManager = project.layoutManager()
for child in root.children():
    if isinstance(child, QgsLayerTreeGroup):
        #print ("- group: " + child.name())
        for subChild in child.children():
            if isinstance (subChild, QgsLayerTreeGroup):
                print ("- group: " + child.name() + "-subgroup: " + subChild.name())
                #turn on subgroup
                child.setItemVisibilityChecked(True)
                subChild.setItemVisibilityChecked(True)
                #open layout
                l_out = QgsProject.instance().layoutManager().layoutByName(subChild.name())
                iface.openLayoutDesigner(layout=l_out)
                mapName= (subChild.name()+".pdf")
                print (mapName)
                base_path = (r"F:\Documents\PDF")
                print (base_path)
                pdf_path = os.path.join(base_path,mapName)
                print (pdf_path)
                #export layout
                exporter = QgsLayoutExporter(l_out)
                exporter.exportToPdf(pdf_path, QgsLayoutExporter.PdfExportSettings())
                #close Layout
                ???

                #turn off layer
                subChild.setItemVisibilityChecked(False)


texts: QGIS performance wtih many unique value widgets in the attributes form

texts:
I'm creating a QGIS 3.8 workflow and made about 15 attributes as unique values generated from the data already in the form but now the info view and attribute table that's used to open almost instantaneously now takes ~30s or more to open and about ~15s when data is updated.

When the options are known is it better to use the Value Map option -this takes longer to setup so is not ideal? or is there some other tweek?


texts: How to make labels' size relative to the map in QGIS1.8?

texts: I want to make my labels' size changed automatically with the scale of the map,so in the menu of "Layer labeling settings",I set the Text size" in map Unit".However,the labels in the selected layer disappeared.I turned off the "Scale-based visibility" and tried to set the size from 0.05 to 2000,but the labels still refused to meet me :P What should  I do?


texts: Error loading AGSMap: "Webmap version is not supported"

texts:
I'm trying to load a map in an iOS app from an on-premise instance of ArcServer.

func initMap() {
    self.mapView.map = AGSMap.init(url: URL(string: self.mapLayerUrl)!)

    self.mapView.map?.load(completion: { (error) in
        guard error == nil else {
            print("Error loading map: ", error?.localizedDescription ?? "???")
            return
        }
        //Success branch never reached
    })
}

This always triggers the error branch and logs the message "webmap version is not supported". This article indicated the map is out of date and needs to be saved. I navigated to the map through the ArcGIS Server site, selected "ArcGIS Online Map Viewer" and saved it. However, that did not solve the issue.

Can someone suggest another way to troubleshoot this?


texts:
# Create a dataframe from sorted JSON dictionary containing desired data
df = pd.DataFrame(data, columns=['lat', 'lon', 'id'])

# coerce dtype 'O' object values to np array native float64 type
df['id'] = pd.to_numeric(df['id'], errors='coerce')
df['lat'] = pd.to_numeric(df['lat'], errors='coerce')
df['lon'] = pd.to_numeric(df['lon'], errors='coerce')

# remove  all rows containing NaNs (no data2 information)
df = df[~np.isnan(osmdataframe).any(axis=1)]

# Create a sorted np array
sorted_osmNParray = np.array(osmdataframe.to_records())

# reconvert to np array so it can be passed to NumpyArrayToTable method
osmNParray = np.array(sorted_osmNParray, np.dtype([('index', '<i8'), ('lat', '<f8'),  ('lon', '<f8'), ('id', '<i8')]))

# Convert NP Array to table
arcpy.da.NumPyArrayToTable(osmNParray, output)

My goal for this script is to extract some information from an Overpass API query and ultimately join it with the table of another feature class I have created.

So far I have completed the following steps relevant to my problem:

Extract the desired data from the JSON dictionary

Put data into a pandas dataframe
convert said dataframe to a sorted numpy array, then back to nparray
export the data in my numpy array to an ArcGIS table via NumPyArrayToTable()

A look at the contents (view from dataframe without nans)

             lat        lon          id
0      48.132752  11.566575  1576559984
1      48.138378  11.576535  4785090864
...          ...        ...         ...
23703  48.166622  11.587858  1601228494
23704  48.166629  11.587299  1601228506
23705  48.166637  11.587443  1601228508
23706  48.166638  11.587735  1601228509
23707  48.166641  11.587711  1601228510
23708  48.166649  11.587303  1601228511
23709  48.166994  11.587971  1601228515
...          ...        ...         ...

A look at the data types of the numpy array

print(sorted_osmNParray_way.dtype)
...
(numpy.record, [(u'index', '<i8'), (u'lat', '<f8'), (u'lon', '<f8'), (u'id', '<i8')])

Finally the contents of the sorted numpy array that is passed to the method

[20532 rows x 3 columns]
(numpy.record, [(u'index', '<i8'), (u'lat', '<f8'), (u'lon', '<f8'), (u'id', '<i8')])
array([(    0, 48.1327519, 11.566575 , 1576559984),
       (    1, 48.1383783, 11.5765354, 4785090864),
       ( 3199, 48.1196275, 11.5380289, 1079422305), ...,
       (23726, 48.1658095, 11.5900809, 3955898986),
       (23727, 48.1626348, 11.5999409, 1324415281),
       (23728, 48.1627603, 11.6000468, 1324415506)],
      dtype=(numpy.record, [(u'index', '<i8'), (u'lat', '<f8'), (u'lon', '<f8'), (u'id', '<i8')]))

This unfortunately returns the following error:

NumPyArrayToTable() RuntimeError: The value type is incompatible with the field type. [id]

This is rather interesting as I have used this method with dozens of other datasets without any problem !
Could someone please tell me what exactly isn't properly working out with the use of 'id' here?


texts: Offsetting Labels to the edge of the map (US MAP) in ArcGIS 10.4

texts:
I have been doing a tutorial in ArcGIS and I'm trying to label cities in California but I can't seem to locate the right setting to bring all the city labels to the edge of the study area. I am using the offset feature in Maplex but I can't get the result I'm looking for. The example in the tutorial shows all cities on the edge of the study area. I'm thinking perhaps this was done manually, but I'm not sure.

Or would there be another alternative to labelling in a somewhat messy map.


texts: Access GeoNetwork thumbnail images on any local server? Relative path issue?

texts:
I've installed GeoNetwork via Docker on an Ubuntu server for use internally with a volume to the host /host/geonetwork-docker:/var/lib/geonetwork_data

I've got 1560 records to import via XML and have thumbnail images for all of them, both on the Ubuntu server GeoNetwork and on a separate windows server used internally as a data repository.

I've tried putting the images in /host/geonetwork-docker/data/Thumbnails/Geo_125k/XXXXX.jpg and using the below path in the XML

         <gmd:graphicOverview>
            <gmd:MD_BrowseGraphic>
               <gmd:fileName>
                  <gco:CharacterString>/data/Thumbnails/Geo_125k/XXXXX.jpg</gco:CharacterString>
               </gmd:fileName>
               <gmd:fileDescription>
                  <gco:CharacterString>Thumbnail</gco:CharacterString>
               </gmd:fileDescription>
               <gmd:fileType>
                  <gco:CharacterString>JPEG</gco:CharacterString>
               </gmd:fileType>
            </gmd:MD_BrowseGraphic>
         </gmd:graphicOverview>

This didn't work so also tried the full path from the host and from the container and still no luck.

I've also tried using the path to the files local windows server http://gdcs1/Metadata/Thumbnails/Geo_125k/XXXXX.jpg which also doesn't work.

Anyone have any suggestions/ideas? Is this a docker issue?


texts: Use of GIS in fixed telecommunication network management

texts:
I want to get an idea how GIS can be used in fixed line telecommunication management?
there are open source GIS tools as well as Proprietary tools as well.
what is the best tool can be used in above ?
thanks.

Chandima


texts:
I want to disable an OpenLayer map's draw pan during activation of a DrawFeature Control. What I am trying to do so far, is overwriting the dragstart method of OpenLayers.Handler.Drag.prototype in the control's activate method and readjusting it in its deactivate method.

However, I am not happy with this since I would prefer to tinker around with the Drag Handler object directly - instead of its prototype - and besides, I am pretty sure there is a more elegant solution which I was not able to find so far.


texts: ArcGIS: Create polygon using file containing point coordinates

texts:
I have a text file containing the projected coordinates of the points that define a single polygon. The X, Y coordinates are separated by Tab

Is there a tool in ArcGIS similar to "Add X,Y Data" in ArcMap, that can be used to automatically create polygons FC, not point FC?

I use ArcGIS 10.1 for Desktop.


texts:
Within a large dataset that contains information on location and attributes of forest plots, I need to identify treatment and reference plots for my research.

I ran into problems with reference plot selection. I'm using a pairwise approach, so each treatment plot is to be paired with a reference plot. The criteria for reference plots are:

 1. of the same forest type as a treatment plot;
 2. located at a minimum distance from a treatment plot.

This means that for each treatment plot, I need to create a subset of plots of the same forest type (a field in the dataset) and then select a plot that is nearest to the treatment plot. Also, I need to make sure that a certain reference plot represented only once.

I have over 500 treatment plots.

How can I do this in in ArcMap?


texts: Are tesselation patterns an effective way of covering a shape?

texts:
I am using Facebook place search to explore business pages on Facebook. If I input a longitude and latitude into facebook's GraphAPI, I can search within a given radius for facebook business pages.

Facebook searches within a circle. Initially, I thought I would be able to get all the businesses within a given circle, so if I made the radius big enough (e.g. the size of New York City), I could get all the businesses (e.g. all the hair salons). The problem is, Facebook limits the number I can return to 100 businesses per search. So even if I am doing a large area, it will cap it off. Therefore, a single circle will not work.

My initial idea was that, for any given shape, one should be able to cover that shape with circles in a tessellation pattern. Here's an example of what I'm envisioning:

Now suppose I have a shape like this:

I assume people have confronted this issue before.

Would this methodology be effective to cover it in circles?


texts:
I have two shapefiles, one is polyline, one is points. I use points to split the polyline. How to record the FID of splited points(two end points each line) in the field of polyline.

First I use split tool

then I execute the following code to update the attribute table of polyline

import arcpy
from arcpy import env
env.overwriteOutput = True
sj = r"in_memory\sj"

inPoints=r"Y:\Desktop\road\city\city_merge1.shp"
inPL = r"Y:\Desktop\road\highway.shp"
arcpy.MakeFeatureLayer_management(inPoints, "city_points")
arcpy.MakeFeatureLayer_management(inPL, "road")

for i in range(0,500):
    expr = " [FID] =" + str(i)

    arcpy.SelectLayerByAttribute_management(inPL, "NEW_SELECTION", expr)
    arcpy.SpatialJoin_analysis(inPL, inPoints,sj,"JOIN_ONE_TO_ONE","KEEP_COMMON")

    ## update the attribute table of road, recording the FID of the city in the attribute table of the road


texts: How to give multiple user input CQL filters to a single layer using Geoserver?

texts: I have a layer in which I need to apply multiple check box filters. For example, I need to filter a layer based on state, city and some other attribute based on user input/ CQL filter. Basically, I need a dynamic filtering based on multiple attributes on the layer using CQL. could any one suggest me an idea on how to do it?


texts:
I'm trying to calculate the area of some census tracts in qgis
I go to vector -> geometry tools -> export/add geometry columns

in the resulting dialog, I calculate using project CRS. What comes out in the attribute table is something like 0.000018 or 0.0000152. I don't know what units these are in; most other information on the internet seems to say that this issue is a map projection issue.

For the record, the current map project is NAD83 EPSG 4269

I do not know what unit this is in. If I knew the unit, it would be pretty trivial to add a new 'calculating' column to the attribute table to determine the area in language I can understand.


texts: How to declare ArcSDE workspace for Geoprocessing tool in ArcGIS 10.1 using ArcObjects and c#?

texts:
What is the proper way to declare an ArcSDE workspace for Geoprocessing tools in c# .net?

I am using ArcGIS 10.1 SDE (MS SQL) on SQL Server 2008. I want to write code so that every night the server will fire a program and compress one of my ArcSDE databases. I thought I could do this in c# using the Compress tool in the data management namespace but have run into a snag.

I am told my workspace is invalid. Esri help docs are not helpful, as they just show generic help and my internet search has lead nowhere.

Can anyone tell me how to define an ArcSDE workspace in 10.1 in net?

Right now I am using @"Database Connections\CME.sde". Here is the error:

Executing: Compress "Database Connections\CME.sde"
Start Time: Fri Oct 12 13:40:47 2012
Failed to execute. Parameters are not valid.
ERROR 000837: The workspace is not the correct workspace type.
Failed to execute (Compress).

What should I type for the in_workspace?

private void button1_Click(object sender, EventArgs e)
        {

            Compress Compressor = new Compress();
            object in_workspace = @"Database Connections\CME.sde";

            Compressor.in_workspace = in_workspace;

            Geoprocessor GP = new Geoprocessor();
            RunTool(GP, Compressor, null);

        }
        private static void RunTool(Geoprocessor geoprocessor, IGPProcess process, ITrackCancel TC)
        {

            // Set the overwrite output option to true
            geoprocessor.OverwriteOutput = true;

            // Execute the tool
            try
            {
                geoprocessor.Execute(process, null);
                ReturnMessages(geoprocessor);

            }
            catch (Exception err)
            {
                MessageBox.Show(err.Message);
                ReturnMessages(geoprocessor);
            }
        }
        private static void ReturnMessages(Geoprocessor gp)
        {
            if (gp.MessageCount > 0)
            {
                for (int Count = 0; Count <= gp.MessageCount - 1; Count++)
                {
                  Console.WriteLine(gp.GetMessage(Count));

                }
            }

        }
    }


texts: Does RapidEye level 3A imagery need any kind of preprocessing before applying classification?

texts:
I just got some rapideye images and I would like to know if there is need of any preprocessing before applying classification to it?
I know that there isn't any automated way of, for example, applying DOS, but my question is: is it actually needed?
Based on the comment below:
I need to apply either min distance or max likelihood, in order to identify built up areas in multiple images (28 to be exact)


texts:
I closed a side panel I had, where I had a list of layers and which I used to switch on and off layers, and I cannot find how to get it back. I have tried all the buttons I could think of, but I cannot see it again.

Can someone help?


texts: Step by step contour creation algorithm?

texts:
I do seismic processing, and recently have been moving into doing contours.

I have watered down my problem so that I can understand from 'first principles'. I have checked out Creating contours for specific elevations as well, since it is similar.

Input: I have a uniformly sampled 2-dimensional (x,y) spatial grid, of size M by N. I have P random points within the grid, making up say, 1% of all the M*N points. Each of those points has a positive z value. All the rest of the points are at z=0.

Output: Based on this information, I would like to know, from an algorithmic / signal processing perspective, the steps needed to make a contour plot.

My problem is somewhat nested:

My questions are the following:

Do I need to first perform a 2D interpolation of my data, so as to fill up all those z=0 values, as a necessary pre-processing step before being able to make a basic contour plot?

1a. If that is the case, then how does one do this 2D interpolation? I am familiar with linear regression, but it seems from my trials that this was not adequate. Nearest neighbours interpolation also does not seem like it would yield good looking contours.
If a 2D interpolation is not a necessary pre-processing step before being able to make a contour, then how exactly, from an algorithmic perspective, are those z=0 values being filled? It seems to me that this is totally necessary first, and if so, I would like to know how this is done.

I am very new to contouring.


texts: Debugging ImportError: No module named arcpy

texts:
I have written a script and am trying to run it in Python IDLE 2.7.6.
I first try to import arcpy, but I get an error message saying:

ImportError: No module named arcpy

Why does it not work?
I am using ArcGIS Desktop 10.2.


texts: Adding new field to existing ArcGIS Report?

texts:
I spent a couple of hours creating report from scratch using report wizard adding only certain fields of feature class that I need to the report (.rlf) window. I save and close report.

A few weeks later I realized that I needed an additional field added to the existing report.

I was only able to achieve this by re-creating report from scratch and adding the missing field at the beginning.

Is there a way to add field to existing report?


texts:
I have a table of points (schools t1) and a table of small polygons (buildings t2). They are within large polygons (p) that represent districts. I am trying to get the average distance from a school to all of the buildings withing the same district. I am really new at this stuff and found a post that seemed to have the answer but I can't quite get it. What am I doing wrong?

WITH SCHOOLS AS (

  SELECT  p.gid, p.cell, t1.name, code2, pub_pri2, enrollment, t1.geom as wkb_geometry
    FROM mesa.schools2 t1
    LEFT JOIN mesa.schldist_cellpop p ON ST_CONTAINS (p.geom, t1.geom)
    ORDER BY p.gid),

DISTANCES AS (
  SELECT p.cell, t1.geom <-> t2.geom as dist FROM
   (SELECT distinct on(t1.gid) t1.gid, p.cell, t1.geom,
      (SELECT t2.geom FROM mesa.school_pop t2
           ORDER BY t2.geom<->B.wkb_geometry
           LIMIT 1) AS other_geometry
       FROM mesa.schldist_cellpop p, SCHOOLS B, mesa.school_pop t2, mesa.schools2 t1)   AS FOO)

SELECT distinct p.cell,
     count(p.cell) over(partition by p.cell),
     round(avg(distance) over(partition by p.cell)::NUMERIC, 1)
     FROM DISTANCES
     WHERE DISTANCES IS NOT NULL

The error I am getting is
ERROR:  missing FROM-clause entry for table "p"
LINE 9:      SELECT p.cell, t1.geom <-> t2.geom as dist FROM


texts: PDAL is mixing Lat and Lon when using tindex

texts:
I am trying to extract polygon representing the extent of my .LAS file. I am using tindex to do so.

When importing the resulting file, it looks like longitude coordinates are set as latitude and vice versa.

This is my code:

pdal tindex create --tindex output.shp  --filters.hexbin.edge_size=1 --filters.hexbin.threshold=4 input.las

Can someone help me understand my mistake?


texts: Catch WMS empty answers with OpenLayers 3

texts:
I was wondering if there is a way with OpenLayers 3 to catch problems with WMS layers, such as the service not answering, or not sending any images/tiles. So that in case of errors with a distant WMS layer (such as OSM WMS or anything), we can switch to another layer.

It is easy to do that by "testing" the WMS GetMap and analyzing an random answer (on a selected BBOX for instance), and analyzing the answer with PHP for instance. But is there a way with JavaScript / OpenLayers ?


texts: Finding variables of QgsPolygon and QgsPolyline?

texts:
now I want to use QgsPolygon (and Qgspolyline,too) class in python

I have a example by using QgsPoint but I don't know how to use in QgsPolygon

following is my Qgspoint example

    feature = QgsFeature()
    qpoint = QgsPoint( float(x), float(y) ) #get qpoint's attitube here
    feature.setGeometry(QgsGeometry.fromPoint(qpoint))#set point's geometry here
    if numpy.isnan(val):
        pass #no assign value means assign NULL value
    else:
        attr_val = round(val,4)
        feature.addAttribute(0,QVariant(attr_val))
    feature_list.append(feature)#append value here

I want to change QgsPolygon(or QgsPolyline) instead of QgsPoint
but I couldn't find the variables of QgsPolygon(or QgsPolyline)

where to find the variables of QgsPolygon and QgsPolyline?

I have another problem here ,

  for s in dataset:
        x=s[0]
        y=s[1]
        feature = QgsFeature()
        qpoint = QgsPoint( float(x) , float(y) )
        polygon=[]
        polygon.append(qpoint)

        feature.setGeometry(QgsGeometry.fromPolygon([polygon]))

but it's still get

        TypeError: QgsFeature.setGeometry(): argument 1 has unexpeced type 'NoneType'

how to fix it?

Have I mistake the meaning of your answer?


texts: Increasing Clickable Area of OpenLayers Polylines

texts:
I'm currently working on a web application where there are OpenLayers geometries that are clickable.  However, since the size of the geometries can be rather small, it can be hard to click on them.  Is it possible to add padding to increase the clickable area without increasing the size of the polyline?

UPDATE: I've since figured out that the way the application works, when a point is clicked, OpenLayers takes the coordinate and tries to find what is at that coordinate.  So I have to find a way to expand the coordinate and send it a few other coordinates to see if it will select what I want.


texts: Why is ArcMap editor deleting some added vertices and simplifying polygon?

texts:
I am using ArcMap 10.0. I am trying to edit (and add) vertices in a polygon.

I start editing, (using "Edit Vertices" tool), and add many vertices in the polygon.

Then when I am done, I save the changes. But, the image gets simplified and I lose many of the additions!

Why is the editor deleting some vertices and simplifying the polygon?


texts: Applying OSM z-order after osm2pgsql with qgis

texts:
I've trying to find a solution to applying the OpenStreetMap z-order in QGIS after using osm2pgsql. Everything is working perfectly, until I add the each layer to the canvas, where the z-order is not being applied.

I've tried searching for solution for the last couple of hours and have yet to find anything that is working.

Here are the steps I've taken thus far:

osm2pgsql -d OSM_Test -l -u  -U 'username' -W -H 'host' -P 'port' -k
'filepath'/oregon-latest.osm.pbf

Note: '<...>' is being used as a placeholder.

Each layer (line, point, polygon, road) is added with a land layer
from Natural Earth.
Rule-based styles are added for each OpenStreetMap road type.

I know that I can combine each road style with the z-order, however, this can create 60 different options for each road type/style.

Where am I going wrong? or am I missing something with a PostGIS function than can applied to the z-order?

I have confirmed that the z-order column is present for each.


texts: Measuring distance from point to closest point feature in GDAL/OGR in same point shapefile?

texts:
I am tying to trying to get the closest point of every point in Python using GDAL module. But some how I did some mistakes.

Here is my sample code:

from osgeo import ogr
from osgeo import gdal

shpfile=r"P:\materials_python\All_GPS_apr17_may5_2016_Corrected_final\All_GPS_apr17Tomay5_2016_Corrected_final.shp"
driver=ogr.GetDriverByName('ESRI Shapefile')
PointsDs=driver.Open(shpfile)
layer=PointsDs.GetLayer()
counts=layer.GetFeatureCount()# get featurecount in gdal
print counts

geomAll=[]
for feature in layer:
    geom=feature.GetGeometryRef()
    print geom.ExportToWkt()
    geomAll.append(geom)
print geomAll[0]
for i in range(0,len(geomAll)):
    dist=min([x.Distance(geomAll[i+1])for x in geomAll])
    print geomAll[i].ExportToWkt()

But it is not showing any result after this line:

print geom.ExportToWkt()

I did not get also any error Message. I do not know what I did wrong.
Does anybody have any idea what is wrong?


texts: lidR - add vegetation classification to unique trees

texts:
I have done a individual tree segmentation using lidR and the resulting LAS file now has been classified so that all points belonging to a tree have a unique treeID. So far so good!
But before writing the LAS file I would like to add High, Medium and Low vegetation to the trees as well. So for example if a tree have a height of let say 5 meter it should be classified as Low vegetation (classification = 3) and between 5 to 10 meters should be classified to Medium vegetation (classification = 4) and so on. I cant figure out how to do it, is this possible using lidR?


texts:
Anyone knows where I can find a complete* database with the historic flood occurrence in Europe, at least since the 70's ?

I am looking for some information on the date, number of floods*, amount precipitation*, as well as the location*** (the region and/or the affected city will be nice).

*complete in time(from the 70's  or 80's until now)
*both informations are not extremely important
**a name will be enough

I am not looking for a database with the geographic coordinates of the flood extent, just a reference for the name of the region/city


texts: Organising attribute table in regard to categories in QGIS

texts:
I'm using QGIS 3.10.2 A Coruna.

I translated a DWG/PDF Plan into GIS where I harmonized every object in regard to topology and correct geometry. The map represents a development plan of a project we are working on. In this development plan we have different categories: buildings, streets, park area/open space, estates.

At first, I had a seperate layer for each categorie, until I realized that everytime a change is made in the plan, It is important to check all the other layers and make sure, that topology is still intact. I found this to be a hassle. So, I decided to put all categories in one layer, that way any changes that occur, will only effect one and one layer only.

However, I have a lot of attributes that need to be implemented in the layer, the problem ist, that some attributes only apply to some categories and not others, and some attributes apply to all categories. For example: The attribute "Name" applies to all categories, but the attribute "use" is only for buildings, where I have predifined values.

My attribute table to growing rapidly and my capacity of keeping track of which attribute is set for which category is decreasing. For my collegues I often have to export only one category, and if I always manually delete unecessary attributes, mistakes are bound to happen.

Is there a possibility in QGIS to only show/export the relevant attribute without me having to do so manually every time I want to look at my data?

I have also considered trying out a database-solution like PostGIS, but I'm new to this, so I'm not really sure how I would organise it there either.


texts: Distance between two latitudes and two longitudes

texts:
I would like to know what is the distance in km between two longitudes and two latitudes.
For example, at the same latitude, what is the distance, in km, between longitude x and longitude x+1 (where x is the longitude in decimal).
Same question for latitude: at the same longitude what is the distance, in km, between latitude x and latitude x+1 ? (where x is the latitude in decimal)
Edit to avoid confusion:
I do not want to calculate the distance between two points with latitude and longitude, I want to know how many kilometers is a difference of 1 for two longitudes at the same latitude and two latitudes at the same longitudes.
For example, distance between longitude 43.0 and 44.0 is XXX km.


texts: Clarification of the "25D" option in ogr2ogr

texts:
[EDIT: In my original post, I had received an error that was based on a typo, so I scrapped that part of the question, but left the questions that still apply.]

In the ogr2ogr docs it says:

-nlt type:
Define the geometry type for the created layer. One of NONE, GEOMETRY, POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION, MULTIPOINT, MULTIPOLYGON or MULTILINESTRING. Add "25D" to the name to get 2.5D versions.

Does this part: Add "25D" to the name to get 2.5D versions apply to both loading things into postgres as well as exporting data from postgres to .shp files?

Additionally, I am assuming 25D to meant that one can have a z value that corresponds to each pair of xy coordinates (as is the case with PolygonZ shapefile types), but that these coordinates may not overlap. Is this correct? What is the intended distinction in this case between 2.5D and 3D?

Thanks


texts: Where is composition tab in QGIS Composer to change page size?

texts: I am trying to change the page size in the QGIS 3.0 Composer but I am not able to see the Composition tab.  I have right clicked an empty space in the toolbar but did not see anything about composition.  I am using Windows 7.  I tried looking at the QGIS docs but they are still in testing.


texts: QGIS and 'Layer to labelled layer' plug-in

texts:
I'm giving a try to this plug-in and I noticed this:
When you use this plug-in on a layer, it adds a lot of field in the layer table.
The field LblColor has a 7 character field, but it needs a 'r,g,b' number that could use more than 7 characters.
In my case I need to write 230,179,251 that are 11 characters...

Is there a way to change this property in the table from QGIS?

If not, I should delete this field and create a new one with a lenght of atleast 12 carachters.


texts:
I want to calculate the Hausdorff distance(https://en.wikipedia.org/wiki/Hausdorff_distance) between two poly lines and i found https://github.com/anitagraser/QGIS-Processing-tools/blob/master/1.1/scripts/hausdorff_distance_pairwise.py ...but i don know how to code in this environment and so i don know how to use it ? I mean in which lines i should insert my inputs and how to write their paths or etc...I would be thankful if you help me...

i will also be thankful if someone recommend me a python code for Hausdorff distance in ArcGIS.


texts: MySQL Equivalent of SQL Server ShortestLineTo()

texts:
In SQL Server I have done the following to find the nearest point of a polygon to a given point outside that polygon:
MyPoint.ShortestLineTo(MyPolygon).STEndPoint()

I'm trying to port this to MySQL (5.7 on Aurora Serverless), but I can't seem to find an equivalent to ShortestLineTo().
Is there an idiomatic way to accomplish this kind of calculation in MySQL?
I searched SE and didn't find much that was specific to MySQL along these lines.


texts: I am trying to generate random points for a survey, but I want the points to fall based on 2 polygons. I have created a polygon fishnet, and have a layer of forested habitats.  I am trying to generate points so that there is one point per grid, and so that the point falls within a forested habitat in the grid.  Is there a straightforward way to do this?


texts: How to rotate labels by the same angle as the map in the composer?

texts:
Does someone know if it is possible to make the labels rotate the same angle as the map in the composer?
I'm trying using the expression attribute( $map, rotation) in the data defining placement, but it doesn't work.


texts: Fetching data from ArcGIS attribute tables using ArcPy?

texts:
I want to fetch/extract data from ArcGIS attribute tables (to set some temporal series) using Python. Because of the structure of my tables and their disposition, I thought of using the "division"

Select ...
from...
Where Not Exists (
   Select...
   From...
   Where...
   And...)
Group By...

I have then to declare or to set cursors for at least two tables and to add the adequate command lines in my script.
I tried this (it doesn't work, I just tried it naively):

import arcpy
rows_n = arcpy.InsertCursor("ndvi_date")
rows_d = arcpy.InsertCursor("Date")
rows_ep = arcpy.InsertCursor("Maroc_Aqua_NDVI_2005_10_decade1_maxPr"
rows_n = rows.newRow()
rows_n.setValue("Date", "Datetime")
rows_n.setValue("NDVI", "Mean")
From nd = rows_n, rows_ep:

 Select "Mean"

 Where Not Exists:

     Select "Datetime"
     from rows_d
     Where Not Exists:
         Select "Datetime"
         from n2 = rows_n
         Where nd.mean = nd2.mean
         And nd2.Datetime = rows_d.Datetime
Group By id_Province

Bellow is the structure of my tables:

So to sum up, what I want is to extract data from my tables, to use  the information in my future scripts


texts: Limit what is shown in the layer panel window under a categorized layer to that which is showing in the canvas extents in QGIS

texts:
I have imported a csv file via 'Open data source manager' > 'Delimited text' which has multiple locations all specified with Eastings and Northings, these are plotted fine and in the layer styling panel I've categorized on a particular header from the csv template and styled accordingly. All works great, however I wanted to know if it was possible and how you would go about it if when you're zoomed in and say only 3 of the 12 locations within this layer are showing on the canvas/screen then only those three locations are displayed in the layer panel?
I know this can be achieved in the print manager via a checkbox setting but I wanted to apply the same functionality to the main 'working' environment too if it was possible...

I'd only want the items circled in red showing when only these sites are in the canvas extent.
Bonus points for anyone who can tell me how I get the icons to be much smaller/not symbol in the layer list without affecting the object on the map canvas!


texts:
Is there some sort of plotting system which lets you do this? (in PHP):
ex. plot(lat,long,data_show_onhover)
My CSV file that I want to plot on the map looks like this:
something,something,something,something,something,Lat,Long,something,something,something,something,something,something,something
So when I hover on the area of lat,long it shows rest of data in array.
I'm using repl.it so there's no root and I want to avoid using google maps because it's paid.


texts: SRIDs do not match for the two geometries in Oracle Spatial

texts:
I am trying to resolve where some points from Point Geometry lying outside the Polygon geometry. I want to filter those points which are lying outside the polygon geometry.

I am getting error where SRIDs doesn't match in Oracle Spatial SQL Query. I have two tables viz, Table 1(POS) and Table 2(Territory). The POS table contains the points geometry information and the Territory table contains the polygon geometry information.

Oracle SQL Query:

select sdo_geom.sdo_intersection(x.geoloc,null,
 y.geom1 ,null ,0.005) from

( select geoloc from Territory b  where  b.territory = 'England') x,
(

select sdo_aggr_convexhull  (
         sdoaggrtype (
           sdo_geometry(2001, 8192, sdo_point_type(p.longitude,p.latitude,  NULL), NULL, NULL),
           0.5
         )
       ).sdo_srid  as geom1
from POS p where p.POS_NAME in(
'LONDON-CENTRAL','LONDON-NORTH','LONDON-SOUTH','LONDON-EAST','LONDON-WEST'
)) y ;

Error while executing the above query:

ORA-13037: SRIDs do not match for the two geometries
ORA-06512: at "MDSYS.MD", line 1723
ORA-06512: at "MDSYS.MDERR", line 8
ORA-06512: at "MDSYS.SDO_3GL", line 1657
ORA-06512: at "MDSYS.SDO_GEOM", line 1458
ORA-06512: at "MDSYS.SDO_GEOM", line 1527
13037. 00000 -  "SRIDs do not match for the two geometries"
*Cause:    A Spatial operation is invoked with two geometries where
           one geometry has an SRID and the other geometry does not have an SRID.
*Action:   Make sure that the spatial operations are invoked between two
           geometries with compatible SRIDs.

I have checked other options like searchWithin but I have rule out because some points were lying on the Polygon border and they were excluded with SearchWithin.

POS Table description

Name         Null     Type
------------ -------- --------------
POS_NAME              VARCHAR2(146)
LATITUDE              NUMBER
LONGITUDE             NUMBER
SE_CODE_NAME          VARCHAR2(106)
MI_STYLE              VARCHAR2(254)
MI_PRINX     NOT NULL NUMBER(10)
GEOLOC                SDO_GEOMETRY()

and Territory Table

Name      Null     Type
--------- -------- --------------
T_ID               NUMBER
TERRITORY          VARCHAR2(78)
MI_STYLE           VARCHAR2(254)
MI_PRINX  NOT NULL NUMBER(10)
GEOLOC             SDO_GEOMETRY()

Types of SRIDs in these table

select * from user_sdo_geom_metadata where table_name in ('POS','TERRITORY')

TABLE_NAME  COLUMN_NAME   DIMINFO                                                               SRID
----------  -----------   -------------------------------------------------------------------   -----
POS         GEOLOC        MDSYS.SDO_DIM_ARRAY([MDSYS.SDO_DIM_ELEMENT],[MDSYS.SDO_DIM_ELEMENT])  8192
TERRITORY   GEOLOC        MDSYS.SDO_DIM_ARRAY([MDSYS.SDO_DIM_ELEMENT],[MDSYS.SDO_DIM_ELEMENT])  8192

Please help me to filter the points which are outside the Territory Polygon.


texts: How to upload a list of 5-digit zip codes to CartoDB?

texts: I'm playing around with the free version of CartoDB right now. Trying to upload both a CSV and Excel document with a list of 5-digit zip codes. I have received a "there was a problem georeferencing your data" error. Any idea why this might be?


texts: Visualizing resource dependencies and relationships underlying EPSG CRS

texts:
I'll go straight to the point; I remember long ago that I was quite shocked the first time I discovered that EPSG:4326 was kind of an "umbrella" definition on top of many others, which often stay hidden unless you specifically need them (which is almost never the case except if you are a geodesist).
Is there a tool to visually represent and explore all the underlying EPSG resources, their relationships and the structure binding them, which altogether make up any of the most well-known CRS, such as the ESPG:4326?
The result I'm looking for should look like what you see in this picture: https://tomroth.com.au/epsg/Capture-12.png
Some kind of visual equivalent for the great tool projinfo.
It's for pedagogic purposes.


texts:
I create a plugin in QGIS using some algorithms from Processing.
When I try to use one algorithm in two different outputs, the variable is not well called and I do not know why.
First call:
Processing.initialize()

params_ndvi = {
'INPUT': str(self.dlg.lineEdit.text()),
'POLYGONS': vl,
'OUTPUT': 'TEMPORARY_OUTPUT'}
feedback1 = QgsProcessingFeedback()
clip_ndvi = Processing.runAlgorithm("saga:cliprasterwithpolygon", params_ndvi, feedback=feedback1)
out_ndvi = clip_ndvi['OUTPUT']

Second call:
params_class_image = {
'INPUT': str(self.dlg.lineEdit_2.text()),
'POLYGONS': vl,
'OUTPUT': 'TEMPORARY_OUTPUT'}
feedback2 = QgsProcessingFeedback()
clip_classified = Processing.runAlgorithm("saga:cliprasterwithpolygon", params_ndvi, feedback=feedback2)
out_classified = clip_classified['OUTPUT']

In this third call, I want to run with the second call output as input but in the final result, in temporary files, the algorithm gives me the report of the first call. I do not know why because the input is well defined:
params_raster = {
'INPUT': str(out_classified),
'BAND': 1, 'OUTPUT_HTML_FILE': 'TEMPORARY_OUTPUT'}
feedback3 = QgsProcessingFeedback()
raster_class = Processing.runAlgorithm("native:rasterlayeruniquevaluesreport", params_raster, feedback=feedback3)
params = raster_class['OUTPUT']


texts:
I need to change stretch of a raster using properties - Symbology to Minimum-Maximum
like image bellow :

Problem is I have about 172 raster, is there any way to make this automaticaly change?
I have arcinfo license


texts: How can I get the popup value to change when I click on another feature, without clicking on a map area where there are no features first?

texts:
I'm working with openLayers 3, which is new to me, and I've been able to add a popup to my map that comes up when you click on a feature in a vector layer. The trouble is the popup content does not update when you click on another feature; you need to click on an area where there isn't a feature to close it, then click on another feature for the popup content to refresh. I'd like for the popup content to update when you click on a new feature, without having to click on an area where there is not feature first. Is there a way I can do this? Here is the code for the popup that I am using:

map.on('click', function(evt) {
  var feature = map.forEachFeatureAtPixel(evt.pixel,
      function(feature, layer) {
        return feature;
      });
  if (feature) {
    popup.setPosition(evt.coordinate);
    $(element).popover({
      'placement': 'top',
      'html': true,
      'content': 'Parcel No:
'+
      feature.get('parcel')+'
'+
      'Owner:
'+
      feature.get('primowner')
    });
    $(element).popover('show');
  } else {
    $(element).popover('destroy');
  }
});


texts: How to display radius while adding element in ArcObject?

texts:
I want to display the radius of drawing circle when user adds this figure on the map.

User clicks on the icon.
Then user begins drawing circle (I process the MouseDown event)

IGeometry geom = Map.TrackCircle();
When user clicks on the map, I draw this circle on the map (I get radius from TrackCircle and center from mouse coordinates).

I want to display radius on the map or in the bottom tooltip (like ArcMap) on the step 2 (TrackCircle).

I think I must process the event OnAfterDraw of the MapControl, but I do not know how (how can I know that user drawing circle and the radius of this circle?).


texts: Extracting data of two overlapping rasters

texts:
I am trying extract the data of two overlapping set of rasters (one, a stack of 35 rasters, all from the same source and the second an elevation raster) to get a data.frame of the values (mean of the values) of each pixel of all the rasters.

The description of the raster stack is the following:

> stack_pacifico
class       : RasterStack
dimensions  : 997, 709, 706873, 35  (nrow, ncol, ncell, nlayers)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -81.62083, -75.7125, 0.3458336, 8.654167  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : F101992.v//ts.avg_vis, F101993.v//ts.avg_vis, F101994.v//ts.avg_vis, F121994.v//ts.avg_vis, F121995.v//ts.avg_vis, F121996.v//ts.avg_vis, F121997.v//ts.avg_vis, F121998.v//ts.avg_vis, F121999.v//ts.avg_vis, F141997.v//ts.avg_vis, F141998.v//ts.avg_vis, F141999.v//ts.avg_vis, F142000.v//ts.avg_vis, F142001.v//ts.avg_vis, F142002.v//ts.avg_vis, ...
min values  :                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0, ...
max values  :                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63,                    63, ...

And for the elevation raster:

> elevation_pacifico
class       : RasterLayer
dimensions  : 997, 709, 706873  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -81.62083, -75.7125, 0.3458336, 8.654167  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : COL_msk_alt
values      : -16, 5164  (min, max)

It is my first time with raster data and I want to extract the data by grids of 1km2 (or less). I know the resolution of both rasters can be coerced to fit into that area requirement, also both dimensions are equal, so the number of pixels per raster is the same.

My question is, can I only merge all the rasters (the ones in the stack and the elevation raster) and extract the data with the confidence that all the pixels overlap (or are in the same place)? Or do I have to create a master SpatialGrid or SpatialPixel object and then extract the raster data to these objects?

Elevation:

#Download country map and filter by pacific states
colombia_departments <- getData("GADM", download=T, country="CO", level=1)
pacific_littoral <- c(11, 13, 21, 30)
pacific_littoral_map <- colombia_departments[colombia_departments@data$ID_1 %in% pacific_littoral, ]

#Download elevation data and filter it for pacific states
elevation <- getData("alt", co="COL")
elevation_pacifico <- crop(elevation, pacific_littoral_map)
elevation_pacifico <- setExtent(elevation_pacifico, rasters_extent)


texts: GRASSGIS Batch Process Command Not Running

texts:
I'm trying to rescale MODIS LST data by running the following command but it gives this error (Image Attached) Kindly help if the syntax is not right.

Command:

for %i (*MOD11A1.A2000065.LST_Day_1km*.tif) Do r mapcalc "%i_rescale=(%i * 0.02)-273.15"

Console log:


texts: Rotating current map extent around center point in OpenLayers?

texts:
What would be the best way to take current extent and rotate it for 90 degrees around center point and it would work for any projection.

Let's say I have a map with n-layers enabled and landscape view. Meaning I can get extent of it. How can I get extent of a rotated view (90° rotation).
Sketch of rotation:

_

Center is fixed and should stay the same and aspect ration should be same as before.


texts: "Create Geospatial PDF" makes PDF dimensions much larger

texts:
I am using QGIS 3.10.9-A.  When I use Export as PDF... and check the Create Geospatial PDF (GeoPDF) box the size (linear - not number of bytes) of the PDF created is much larger than the desired paper size.  This happens even if the layout is simply a blank sheet of paper.  If I make another PDF without checking the Geospatial PDF option the size is the same as desired paper size.
For example, if I specify the layout to be on Letter-size (11 x 8.5 inches) paper and then open the PDF in Adobe Acrobat, the "Document:" dimensions shown in Acrobat when preparing to print with page scaling set to None are 45.8 x 35.4 inches, 4.16 times the paper size.
I can print the PDF on Letter-size paper using page scaling set to Fit to Printable Area or Shrink to Printable Area (both show a Zoom of 23%) but to print to accurate scale I prefer to set page scaling to None.  Also, if combining PDFs of QGIS layouts with PDFs from, say, MS Word, the reader will have to rescale to read.
Any suggestions?


texts: Automatically update fields in polyline based on points snapped to the beginning and end

texts: I have a bunch of points with unique ID's and polylines joining them. I need to update the fields on the polylines to have the ID of the point snapped to the beginning and end of each line. How can I do this in MapInfo Pro 2019?


texts: Getting KML from blob file of Oracle using ArcPy?

texts:
How do I get a kml which is stored in a blob file of oracle using ArcPy?

i try get data or kml files contains in a database oracle  in a field of a table, using cursors options  of arcpy.

import arcpy, os
from arcpy import env
source =r'C:\connectionFiles\file.sde\tableoracle'
cursor = arcpy.SearchCursor(source)

The conection works correctly,  but the cursors do not work well for this type of information.


texts: Multiplying z values in linestringz shapefile with a constant

texts: I have National Levee Data (NLD) shapefile for the entire US in a linestringZ format. The unit of Z values in the shapefile is in feet, but I need it in meters. So, I want to multiply all z values with 0.3048 to have the elevations in meters. What tool could I use to do this multiplication? Or any other ways to do this conversion.


texts: Aggregating Individual District Assignments to Census Block

texts:
Working on a solution that may work but figured I'd ask the community to see if there's a smarter way to get this task done.

I have a voter registration file with district assignments for each voter. I have a census block assignment for each voter linked on their unique identifier. 95% of those census block assignments are accurate.

Is there an easy way to "build up" to the census block to confirm accuracy? For instance, if I know voter x lives in County A, Precinct 3, and Legislative District Y---how can I programatically ensure that that voters assigned census block falls into the list of possibilities for that combination of districts?

Using this process, I could potentially narrow down the list of possible assignments to 3-5 blocks for the 5% of voters I'm not confident have correct assignments and begin cleaning my data set up better.


texts: How to calculate the ratio between DPI and MB in QGIS?

texts:
I need to make maps that do not exceed a limit of 20 MB.
Up until now I worked with DPI values, that I choose when I export maps from the print composer

Maps with the same DPI have a different  MB of, a fact that I think is link to the complexity of the project. As an example a map with a HD raster will be with higher MB relatively to a similar map without the raster, both maps with the same DPI.
My question is: Can I choose the MB of the export image like the DPI?


texts: Getting Uncaught SyntaxError: Unexpected token { error when using downloaded ArcGIS API in PC without internet)

texts:
I am trying to implement ArcGIS API offline. I have downloaded the zip file from the site and followed the following changes, changed the baseUrl to my Apache IP copied the files to /var/www/html/ and replaced all occurances of https with http in dojo.js and init.js files (My server users http).
After that, I used the HTML code given in the installation isntructions to check if the installation is successful or not.
I tried opening/running the HTML in machines with internet, it worked perfectly, but when I try to run it in machines without internet connect , it shows "Uncaught syntaxError".

I feel like it has something to do with "http://www.w3.org/2000/svg" in line number 68 in init.js ,but i dont know how to solve the error.
Edit 1 :- I have tried to understand the error and even remove the code which is causing the Extension. It appears that using catch like catch{} is causing this error. As far as I know, catch(e){} should have been the correct syntax.


texts: Opening HDF5 satellite images in GDAL with Python?

texts:
I want to read and open HDF5 satellite images by python for batch processing.

How can I read and open both image bands and meta data?

For batch processing I need to geo-reference and read IR bands. Then clip them and calculate Brightness temperature for IR bands based on some formulas.


texts:
We have ESRI ArcMap 10.0 (build 4400) ArcGIS Desktop 10 Service Pack 5. SQL Database 10.0.4000.

We are trying to Compress our Versioned tables but for some reason it's not pulling the changed records into the SQL Tables.

I have checked the SDE_compress_log and can see the Compression was successful but the start_state_count is 5 and the end_state_count is 5 so I think this means nothing was released from the delta tables? Is that true? I also noticed that the end_state_count never seems to be zero. Why is that?

I checked the SDE_versions table and there is only 1 record with the latest state_id, however there are multiple updates that still need to be compressed. There are no related replicas. The compression worked last prior to July 7th 2016 but now nothing is updated. Should the SDE_versions table contain a record for each pending compression or just the last one?

I did some research and created the following query. It returns 5 rows. One for each pending update.

SELECT mv.[state_id],
   mv.[registration_id],
   tr.table_name,
   s.creation_time,
   sl.lineage_name,
   sl.lineage_id,
   lm.time_last_modified
FROM [sde].[SDE_mvtables_modified] mv
JOIN sde.SDE_table_registry tr ON mv.registration_id = tr.registration_id
JOIN sde.SDE_states s ON mv.state_id = s.state_id
JOIN sde.SDE_state_lineages sl ON s.lineage_name = sl.lineage_name
JOIN sde.SDE_lineages_modified lm ON sl.lineage_name = lm.lineage_name
JOIN sde.SDE_versions v ON s.state_id = v.state_id

Can anyone explain what might be happening? It's like there is an unfinished transaction that is preventing the release of the delta data.

This is what I see in my SDE_state_lineages table. Is this normal?

lineage_name lineage_id
0            0
15908        0
41477        0
15908        32512
41477        32512
15908        41276
41477        41276
41477        41548
41477        41552

This is the SDE_versions table

name    owner   version_id  status  state_id    description parent_name parent_owner    parent_version_id   creation_time
DEFAULT sde 1   1   41572   Instance default version.   NULL    NULL    NULL    2011-07-07 14:03:14.000

The table registration_id is 292 so I looked in a292 and there are 2300+ rows in there. Should this "a" table empty out once the compression has successfully pushed the data to the destination SQL table?


texts:
I want to align text to a specific angle in QGIS. Both the text and the angle are specified in a table (column: street name / column: angle) When I generate the labels, at the moment all the labels are placed horizontal. I want the labels (street name) to align to the specified angle from the angle column.

Any ideas on how to do this? Maybe using an expression?


texts: What are the options for publishing maps online, with support for high-res aerial imagery?

texts:
The place I am working wants to publish their maps online - however the only system they have is mapshare, which isn't supported by ESRI anymore. I've suggested DEKHO but I think they are looking at not spending much money. Atm we are going with layered pdfs however I can't export any map with aerial imagery as its 15cm resolution and the pdfs end up blank or only with part of the imagery on it.

Could anyone suggest any other options? I was thinking along the lines of a python application but I am not that good with python and I wouldn't know how to do things server side. The issue I have with layered pdfs is that it doesn't allow people to choose their own extents.


texts: QGIS 2.4.0 doesn't support multiple types of geometry in GML file

texts:
I'm testing GML 3.2.1 data according to INSPIRE Data specification (version 3.0). Data for theme Cadastral Parcels shall have both polygon geometry and definition point. When I download data from my WFS service to QGIS 2.4.0, it renders automatically as points. I'm used to it. In older versions of GIS, all I had to do was overwrite geometryType attribute from 1 to 3 in automatically generated .gfs file. Yes, it's a bit complicated, but it worked.

In QGIS 2.4.0 it doesn't work. After rewriting .gfs file, I got no geometry at all. If I want to get polygons, I have to download GML data and erase all referencePoint elements. So my question is: does anybody have any idea how to visualize GML data with both point and polygon geometry as a polygon in QGIS?

This question may rather be asked to QGIS developers, but at first I'm trying existing solutions.


texts: Semi-Automatic Classification Plugin Syntax Documentation?

texts: I am want to integrate SCP module for DOS corrected into my script. However there appears to be no documentation of its python syntax structure available anywhere?


texts: Order of reading in Row Ordered USGS DEM (.dem) file

texts:
I am making a .dem file parser and so far https://web.archive.org/web/20180612184054/https://nationalmap.gov/standards/pdf/2DEM0198.PDF has proven very useful in providing the information needed.
However, on page 12/92 it mentions Row ordering. So far no example I've found has different setup than 1 column per B record, which you read from South to North and from West to East.
I would like to know the order of how to read a differently ordered file (if it actually exists).


texts: Satellite signature of reindeer carcasses

texts:
In a recent incident (august 2016) in the Norwegian mountains 322 reindeer (rangifer tarandus) were killed by lightning. The carcasses position are documented, although only by point data and not polygons. No vertical photos exist as yet. The carcasses will be left to decay. Heads were removed to do studies on the potential presence of chronic wasting disease in Norway.

Given that we get good vertical imagery for the actual area and some of the surrounding area, I imagine we will have a good point of reference to use the incidence to learn about other similar incidences should they occur. The rangers might not always have time, resources or be lucky enough to find such incidences.

I imagine that indicators of such instances could be change in vegetation due to bird droppings and nutrients in the ground due to carcasses.

I am thinking that researchers might use this incidence for a study where one looks at the possibility to:

find indicators of similar incidences using satellite imagery (by
current signature)
find historical incidences by learning how such an incidence changes the vegetation

Would securing drone photos give us a good point of reference for future studies?


texts: How to create a new attribute indicating which polygon from another layer contains each polygon on a layer

texts:
I have two vector layers composed of polygons. The polygons of the first layer are contained in exactly one polygon of the second layer.

How do I create an attribute on the first layer with the ID of the polygon containing it?

The layers are shapefiles.

I have tried to use the field calculator to create the new attribute, but I could not figure out how refer to the geometries of the second layer.


texts:
I am using map tips to display geo-tagged images on QGIS 2.16.3, but I'm running into problems trying to control the overall size of the map tip pop-up frame.

I have followed the advice on other questions (e.g. Change size of frame for qgis map tip containing text) but it doesn't solve my problem...

This is what I'm using:

[%CONCAT('')%]

[% "DESCR" %], [% "image" %]

This shows the photo, with text below. However, the 'div' doesn't expand beyond about 400px wide so it limits how large the photo can be displayed - the map tip window then shows scroll bars to view the 'overflow' image and text below (as it does if I set the image size larger than the div, as expected) - I haven't found any documentation suggesting there is a limit to the pop-up map tip, but maybe it exists somewhere?

Before anyone asks why I need a large map tip and it would use up lots of my screen, I know that - I have specific desire to see the images as large as I can rather than just a 'thumbnail', and map tips are a great way of introducing an extra dimension to the map.

As a related follow up (which I can post separately if that helps) - is it possible to change the x,y location of the map-tip pop-up. It looks like it appears below-right of the point that I hover - can it be centre-centre, for example. Obviously when I've got a large pop-up, it would be good to be able to see all of it!

I've also found that if I have another layer that uses a similar map tip, but without the 'div' statement within the same project, it displays as (approx) 100px high and 400px wide. Odd behaviour... I know I can get round that by setting a 'div', but that's not the point!


texts: What is X Y coordinate in this data set?

texts:
Anyone knows what the x y coordinate in this crime data set is?
https://data.cityofboston.gov/Public-Safety/Crime-Incident-Reports/7cdf-6fgx

I think it is not lat, long. How do I convert it to lat long?


texts: Add GTFS to Network Dataset; stop_times.txt data is not transfering

texts:
I am using the Add GTFS to Network Dataset toolbox from ArcGIS.com with a publicly available feed.

I successfully imported this feed once using the tools.

however, when I am trying to repeat the process, the stop_times table in the GTFS.sql table is empty, although the fields were loaded.

Any advice on how to correct this behavior?


texts: Polygon features to parts of a multipart feature with chosen attribute source (QGIS)

texts:
In QGIS 3.16.  I have a shapefile (polygon) layer with an attribute table.
What I'd like to do is select certain polygon features and turn them into parts of one multipart polygon feature,
without affecting the other, unselected features on my layer.
If possible, I'd like to do this

in a way that allows me to select which of the multiple selected features to keep the attributes from;
in a way that does not create a new layer each time (I want to do this many times for my layer) but simply collapses the selected features into one line (with the chosen attributes).

If possible, though this is less important, I'd also like to do it:

by physical selection (clicking) rather than grouping by attribute.

I believe that this is possible in ArcGIS, but in QGIS the existing methods of turning 'features' into 'parts of features' that I am aware of, do not allow one or more of these criteria.
For example,

"Collect geometries" automatically picks one feature from which to pull the attributes that it will assigned to the new, multipart feature (and it creates a new layer).
"Aggregate" has a range of options for determining the attributes of the new, multipart feature (collect, sum, concatenate...) but as far as I can see it does not allow simply selecting one of the source features for the totality of its attributes. (And it creates a new layer.)
"Promote to multipart" creates multipart geometry but does not actually turn features into parts of features.

Is there a plugin for this?


texts:
I have an Excel file with about 70 addresses. Is there a way to import these addresses into ArcMap to show up as points in the correct location?
Or do I need to pick a coordinate system?

If so, how do I decide which one?

Basically, I need to find out the distance from each of these addresses to two locations.


texts:
I'm attempting to produce tiles for a tif with the following call to gdal2tiles.py:

gdal2tiles.py --profile=raster C:\Users\kcwil\media\bmng_qsc\0\frontside.tif

and I'm getting the following error:

File "<stdin>", line 1
gdal2tiles.py --profile=raster C:\Users\kcwil\media\bmng_qsc\0\frontside.tif
                               ^

"C:\Users\kcwil\media\bmng_qsc\0\frontside.tif" is a valid file produced by gdalwarp.
As I understand it, gdal2tiles takes its parameters in the form : [options] input_file [output_dir] and what I am doing appears to be similar to examples I've seen.

What am I doing wrong?


texts: Getting Error 1: Failed to create directory for shapefile datastore?

texts:
On page 315 in QGIS training manual (using QGIS 2.18) I am attempting to use the Polygonize from raster to vector tool to vectorize my raster (of forest strands) to polygons for Exercise 14.3.2, I am unable to do this.

An error message keeps popping up when I try and use the tool;

  Error 1: Failed to create directory rautjarvi_green_polygon for
  shapefile datastore.  Traceback (most recent call last): File
  "C:\PROGRA~1\QGIS~1.18\bin\gdal_polygonize.py", line 197, in dst_layer
  = dst_ds.CreateLayer(dst_layername, geom_type=ogr.wkbPolygon, srs = srs ) AttributeError: 'NoneType' object has no attribute 'CreateLayer'

The raster layer has band values of zero where there is no information (no forest boundary). Would this possibly be the cause of the issue, as polygons are only needed for the forest boundary?


texts: Categorising point layer twice by size and form in QGIS

texts:
I have data points from a CSV-file that I opened in QGIS.
The points illustrate different power plants of solar, wind, geothermal, coal, oil, etc with different capacities within each class.
I already categorized the symbols by opening the Layer Properties...> Categorized and selected the value I wanted to categorize. In this case, I first wanted to be able to distinguish which type of power plant it is. Each type has a different color.
However, I would also like to categorize such points by the energy capacity of each plant. For this, I would like to change the size of the points if their capacity falls within a range.
How can this be done in QGIS?


texts: ArcPy: check in and out of Advanced/Standard concurrent license

texts:
Is it possible to check in and out of a concurrent license straight from ArcPy. For example, sometimes I need to use a buffer right or find duplicates. Both tools require ArcPy advanced or Standard. Although we have a large amount of Basic licenses, we have a small amount of concurrent standard/advanced licenses.

Is it possible to do the following (Pseudo code):

# ...

# Pseudo code:
if license == 'BASIC':
   check_in('STANDARD')
   func_that_does_buff_right(data)
   check_out('STANDARD')
   check_in('BASIC')


texts: Near analysis by matching attribute in two different layers

texts:
I need to perform a near analysis to assign cases in one feature class to their nearest match on an attribute in a different feature class (ex: Tutors and Tutees need to be paired by same language spoken). I have not been able to find a way to do this. The closest I have reached to simplifying this process is to perform a split by attribute on the values for the matching variable. So instead of having all tutors in one feature class, I now have n language number of feature classes containing tutors and n language number of feature classes containing tutees. I then plan on performing a near analysis between each feature class. However this is extremely tedious and may result in 100+ near analyses needing to be performed because there are other criteria I will need to match independently on in a similar way.

Is there a better approach?

Much thanks in advance.

Edit: The near must be performed between two different layers to preserve the different roles of individuals being matched. I am new to learning modelbuilder so any ideas on conditionals, or iterators that could determine the "nearest" match on a shared attribute for points housed into two different layers would be much appreciated!


texts: Spatio-temporal kriging - useful if not interpolating to new time points?

texts:
I have data measuring a variable Z say at several locations on a fixed network (each location at a fixed longitude and latitude).  values of Z are measured daily at the same time.

At the moment I am employing ordinary Kriging to interpolate values for each day (ie. run the procedure separately for each day). I am using the gstat package in R.

However, when I look at the autocorrelation function of Z at a given location I see values are correlated up to a lag of about a week. Can I use this additional information to improve my spatial interpolation estimates using spatio-temporal Kriging?  (Just to be clear, I want to interpolate spatially, but to the same times that the data is measured at.) Or is it the case, that by interpolating to the same times as in the data I am essentially using the one dimensional slice of the space-time variogram that corresponds to zero time lag, which means I am just doing ordinary kriging?


texts: How to measure the distance of multiple points to a line

texts:
I have two shapefiles in QGIS to which I am trying to find the distances between. I have a shoreline which is a vector line and then I have hundreds of data points in the ocean.

What I am trying to find is measuring the distance of each point to the corresponding place horizontally on the shoreline thus giving me a distance of each point from the shore. I have tried the distance matrix plugin which seems to only deal with point to point distances and have tried the mmqgis hub distance which only measures distances to a central fixed point on the line. Is there a plugin or way to find the distances I am looking for?