texts: How to determine if a Polygon is inside another Polygon's hole

texts:
Given two valid polygons P and Q, what's the cheapest / fastest way to determine:

If polygon Q is inside one of polygon P's holes (interior rings), or vice versa. (see illustration below)

To the best of my knowledge, the DE-9IM predicates don't directly help here, since the polygons are disjoint (FF2FF1212) according to the matrix, so checks for containment, intersection, crossing, and within aren't useful.
The best I've come up with is:

check whether the polygons intersect

if they don't, construct polygons from their bounding boxes, then check whether P_bbox contains Q_bbox, or vice versa.
One polygon's bounding box will always have to be strictly smaller than the other.

But I'm wondering whether this might catch some other arrangement that I haven't thought of.


texts: Toggling multiple layers with one click? Mapbox GL JS

texts:
How would I go about grouping the two layers in the sample code below so they could be toggled on and off together with one button?

    <meta charset='utf-8' />
    <title></title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.22.1/mapbox-gl.css' rel='stylesheet' />


texts: Google Earth Engine: Apply ee.Reducer to all image values

texts:
Background

I would like to return percentile values for a multi-band image, using ee.Reducer.percentile().

For example, the image could look like this:

var image = ee.ImageCollection('LANDSAT/LT05/C01/T1')
  .filterDate('2008-01-01', '2008-03-01')
  .filter(ee.Filter.eq('WRS_PATH', 44))
  .filter(ee.Filter.eq('WRS_ROW', 34))
  .select('B1', 'B2')
  .toBands();

print(image);

In my case I am stacking an imageCollection (every image is one year), resulting in 18 bands:

var year_stack = imageCollectionYears.toBands();

In the R environment, I would use quantile(values(image), p = seq(0, 1, by = 0.2), na.rm = TRUE), which would return 6 quantile values, taken from all cell values in the raster image, regardless of layering.

Tried so far

I have applied the reducer to the image, but this returns percentile values for every band separately, rather than for all values in the image:

var percentiles = image.reduceRegion({
                  reducer:ee.Reducer.percentile([0, 20, 40, 60, 80, 100]),
                  bestEffort:true
  });

print(percentiles);

In my case, each band is an individual year, and the code returns percentile values for each year.

I have also tried to apply the reducer to my imageCollection, which returns an image with 6 bands.

var test = imageCollectionYears.reduce({
                  reducer:ee.Reducer.percentile([0, 20, 40, 60, 80, 100])

  });

I am struggling to retrieve the numbers of the percentiles in the returned image test, so I don't know if this worked (e.g. print(test.get("b1_p20")) prints null):

Question

How do I apply the reducer to all values, regardless of what band they are in, and return 6 percentile values? My aim is to get overall percentile values over all years, which I would then use for further analysis.


texts: gdal_translate 16bit to 16bit truncating whites

texts:
I have a .CUB map from the USGS, the Mars MOLA image, which is about 46000x23000 pixels and 2GB. I am trying to convert it from .CUB to .TIF but I want to keep it 16 bit.

I ran:

gdalinfo -stats Mars_MGS_MOLA_DEM_mosaic_global_463m.cub

Output:

STATISTICS_MAXIMUM=21241
STATISTICS_MEAN=-720.52571146212
STATISTICS_MINIMUM=-8201
STATISTICS_STDDEV=2976.6805556461

Then I ran:

gdal_translate -scale -8201 21241 1 65536 -a_nodata 0 Mars_MGS_MOLA_DEM_mosaic_global_463m.cub Mars_Bump_16bit.tif

The output image, the four volcanoes on Mars, are not white, they are flat and gray, it's like the tops of the volcanoes were outside the range and were truncated.

If I run:

gdal_translate -scale -8201 21241 1 255 -ot byte -a_nodata 0 Mars_MGS_MOLA_DEM_mosaic_global_463m.cub Mars_Bump_8bit.tif

It generates an image where the volcanoes are NOT truncated.

What am I doing wrong for 16bit to 16bit?


texts:
I am using pgr_dijkstra for finding the shortest path between two given nodes. The network map is actually a grid based one and I have created all the necessary columns and working fine. But I do not understand how does the id2 (which is edge) given?

The values given for id2 (e.g., 59721224, 59721496 ...) can not be found any of my tables. So how does this id2 (edge) values are given? Every single time I query pgr_dijkstra, id2 values are changing but none of the other columns are changing.

This is how my wp_norway_network1 table looks:


texts: CSV file does not display points in QGIS 3.6

texts: I insert a CSV file to QGIS 3.6, but something is wrong with the CRS. Firstly, I set the CRS of the project that I have ( USER 100000 * Παραγόμενο ΣΑΣ (+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396000 +b=3396000 +units=m +no_defs), but the points are not on the CTX Mars mosaic I want. When I zoom in to layer the points are in another location. When I change the CRS to GCS_Mars_2000 or WGS 84 and I zoom to layer again, I am transferred to the correct CTX location, but I do not see any points anywhere. Even when I zoom out to see if the points are somewhere else on the map, they are not.


texts: how transforms coordinates etrs89 in wgs84 with javascript

texts: I'm programming a web page with Visual Studio and I'm trying to make a javascript to use my page. The script has to transform the coordinates X, Y etrs89 WGS84 coordinates that can read Google. I downloaded proj4js-2.4.3 to get the proj4.js and proj4-src.js files and import them into my project. Does anyone know how that script should be written?


texts: Publishing mosaic dataset to image service

texts:
I'm having problems publishing a mosaic dataset to an image service. While I can publish a service, the result does not carry over the pixel values (every cell is valued at 1.134274E+38). Strangely enough, the thumbnail image in arcCatalog shows the correct image, not sure why I can't visualize it any other way.

Some specifics about my mosaic dataset:
it has a time dimension,
contains 12 rasters,
1488 x 802 grid/ single band,
spatial reference: Albers equal area (projected),
pixel type: floating point,
I have ArcGIS Image Extension installed,
I didn't register the underlying data, instead transferring everything (12 rasters) to the server

I've included some screen shots of the process I went through (see below).  Too many steps involved for someone to 'spot' a problem, probably.  I just want to see if someone has had a similar issue and what they did about it. Also, if someone has found a useful guide, I'd love to hear about it. I had used: http://resources.arcgis.com/en/help/main/10.1/index.html#/Publishing_image_services/01540000045v000000/

screen shots:
What I began with in mxd

Making the image service through Arc Catalog

Correct image appears in catalog thumbnail
results not carrying over to preview, services or mxd


texts: OpenLayers and Ionic 3 | tileLoadFunction not working

texts:
I am building a simple map app with the Ionic 3 framework.  I am using OpenLayers 4.0.2 as the map api.

It has been going great until I tried implementing a custom tileLoadFunction for the basemap tile layer used in the app.  In my case, I am attempting to get tiles from an IndexedDB.

 tileLoadFunction: function(tile, src) {
     let imageTile = tile as ol.ImageTile;
     let img = imageTile.getImage() as HTMLImageElement;
     let tileCoord = imageTile.getTileCoord();
     tileCoord[2] = -tileCoord[2] - 1;
     self.mapDataService.getBasemapTile(tileCoord).then(b64 => {
         if (b64) {
             img.setAttribute('src', b64);
         } else {
             img.setAttribute('src', src);
             self.mapDataService.cacheBasemapTile(tileCoord,src);
         }
     });
 }

This is fairly common requirement, and I have seen dozen of implementations online - but not in TypeScript...  I can't just do tile.getImage().src = b64; or else I get complile errors.  So I cast to ol.ImageTile and HTMLImageElement.

I am 100% positive that the images returned by b64 or src variables are valid images. But the basemap tile is always blank???  Is there a better way to do this in TypeScript?


texts: Setting output name to be same as input name in ModelBuilder using Parse Path?

texts: Can you comment on my model as to why the output name is totally different from input name?


texts:
I have implemented MapServer in a Windows environment and I get a map from the following URL

http://localhost/cgi-bin/mapserv.exe?
  mode=browse&
  template=openlayers&
  layers=all&
  map=C:/ms4w/apps/osm/basemaps/osm-google.map&

I need to mark multiple locations on a map, (using map icons) and need to update the map periodically, how do I achieve this?

I guess I have to add a new layer using OpenLayers but can't figure out where I need to make my changes.


texts: Need help: Sorting elements of a combo box in ArcGIS Java script API (ver. 2.8)

texts:
I am using sorting elements in combo box in ArcGIS Java Script API (v.2.8) & for sorting I am using “ fetchProperties="{sort:[{attribute:'name', descending:false}]}".
All entries in combo box are sorting excluding “ALL” entry.

So I am not able to understand that why “ALL” is not sorting ascending. It is displaying at the end but when make “descending: true” then its displaying first.

I want to put “ALL” in the beginning of the combo-box. Therefore, any help will be great.

Please check the screenshot below: enter image description here

Image :

Code is :

  //combo box item
  function populateList(results) {
    //Populate the dropdown list box with unique values
    //alert("sss:  "+results.length);
    var zone;
    var values = [];
    var testVals={};

    //Add option to display all zoning types to the dropdown list
    //values.push({name:"ALL"})

    var features = results.features;
    //values.push({name:'ALL'});
    dojo.forEach (features, function(feature) {
      zone = feature.attributes.TOPIC;

      if (!testVals[zone]) {
        testVals[zone] = true;
        values.push({name:zone});

      }
    });

    var dataItems = {
           identifier: 'name',
           label: 'name',
           items: values
    };
    var store = new dojo.data.ItemFileReadStore({data:dataItems});
    dijit.byId("mySelect").store = store;
    values.push({name:"ALL"})
  }


texts: Downloading an unzipped File Geodatabase from USGS

texts:
I am attempting to download a file geodatabase from the pubs.usgs.gov site.

The file db is DDS69ff.gdb

The browser sees the file as a directory so I am having trouble downloading the complete file as a geodatabase accessible in Arcgis 10.

I've tried to retrieve the file in wget and I've had no luck.


texts: Determine sun elevation without header file?

texts: I have some AWiFS and LISS-III data that I'm trying to convert from at-sensor spectral radiance to top-of-atmosphere reflectance.  The problem is, the data didn't come with any header file.  All I have is an XML document but it doesn't provide any sun-earth geometry.  I could use something like: http://aa.usno.navy.mil/data/docs/AltAz.php knowing the acquisition date, but I can't determine the time of acquisition.  I'm using ERDAS Imagine 2013 and the data-sets were obtained from Bhuvan Gateway to Indian Earth Observation at http://bhuvan-noeda.nrsc.gov.in/download/download/download.php


texts:
Within Google My Maps (https://www.google.com/mymaps) There is an option when exporting the data out called:

"Keep data up to date with network link KML (only usable online)"

I want my Google My Maps to update whenever I update my KML file on Google Drive and I  assume that's what the above does? But why is it on the 'export KML/KMZ' button rather than when you import the data originally?

When you click the 'learn more' button it takes you to a very unhelpful page: https://support.google.com/mymaps/answer/3109452?visit_id=636743436468746992-1430392972&p=lite_kml&rd=1

Kind regards,


texts: Using GRASS plugin tools from OSGEO4W command line

texts:
I am trying to do the interpolation of GRASS vector shapefiles using GRASS plugin tools.
My objective is to run interpolation tool in GRASS plugin from OSGEO4W command line, which will help me automate the process in my project.

I have found an article on how to do this, but it seems to work for QGIS tools and I don't know how to access the GRASS tools in a manner similar to that mentioned in this website. The website URL is as below:

http://nathanw.net/2011/09/27/generating-contours-using-gdal-via-shell-or-qgis/

So if anyone has an idea as to how to do this, please help me out with this.


texts: How can I extract bbox of polygon geometries in a shapefile using geotools

texts:
I have a shapefile which has building footprints as polygon features. I want to read shapefile and extract these bbox of these polygon features using geotools. So far I read shapefile but I don't know how can I access geometry and attributes of a shapefile records using geotools. Here is my code; but I don't know how can I do the rest of work for my purpose.

package org.geotools;

import org.geotools.data.DataStore;
import org.geotools.data.DataStoreFinder;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureSource;

import java.io.File;
import java.io.IOException;
import java.net.URL;
import java.util.HashMap;
import java.util.Map;

public class Tiler {
File file = new File ("C:/Users/Ziya/Desktop/Kasustu/footprints.shp");
Map<String, URL> map = new HashMap<String, URL>();
public Tiler() throws IOException {
    map.put("url", file.toURI().toURL());
}

    DataStore dataStore= DataStoreFinder.getDataStore(map);
    SimpleFeatureSource featureSource =
    dataStore.getFeatureSource(dataStore.getTypeNames()[0]);
    SimpleFeatureCollection collection = featureSource.getFeatures();


texts: Avoiding lock issues when running ogr2ogr concurrently

texts:
I have a Python script that does some stuff and then updates a big geodatabase with the resulting shapefile. I am running this in parallel so occasionally the ogr2ogr command will fail due to lock issue. I am not surprised since multiple shapefiles could be being written to the same gdb at the same time:

ogr2ogr -nln ATL__ZonalStats -a_srs EPSG:4326 -t_srs EPSG:4326 -update -append -f "FileGDB" ATL.gdb ATL__zonalStats.shp

ERROR 1: Failed at writing Row to Table in CreateFeature. (Cannot acquire a lock.)
ERROR 1: Unable to write feature 3128 from layer ATL__zonalStats.
ERROR 1: Terminating translation prematurely after failed translation of layer ATL__zonalStats (use -skipfailures to skip errors)

Is there a way I can get around this issue? Like possibly test whether or not a lock cannot be acquired and wait until it can in order to run the ogr2ogr command? I do not know of a way to do this and cannot find an answer (that is not arcpy-related) on here.

Alternatively I could create the big gdb once all the individual output shapefiles are completed but I would really rather not do that.


texts: Show only some attribute values in a table using PyQGIS

texts:
From Showing attribute table in form using PyQGIS I know how to show a table with all the attributes using this code:

cache = QgsVectorLayerCache(layer, 10000)
model = QgsAttributeTableModel(cache)
model.loadLayer()
table = QTableView()
table.setModel(model)

Now I'm wondering how to show not all the attributes, but only the ones I want or, if it's possible, do the same thing with the classic attribute table.


texts: ESRI Geocoder Service Different Name For Single Line Input

texts:
I have two ArcGIS Geocoding Services but they have different Service Descriptions on Single Line Address Field: part

For This Service I have SingleLine

but for This One it is Single Line Input

Can you please let me know why they are different and How I can configure them on Server or locator side?


texts: Drop Markers and MarkerCluster together from database

texts:
I've been working on my map which gets it's details from a MySQL database.
I use MarkerCluster Plus to cluster my markers but would like to add the drop animation to it.

I found this thread, but cannot get it to work with my database.
Using Drop Markers and MarkerClusterer together?

Would anyone please be able to help me figuring out how to add the drop animation to my code? My test map is here: http://mymap.mooo.com

Really appreciate any help or pointers.

 var customIcons = {
  BOO: {
    icon: 'http://labs.google.com/ridefinder/images/mm_20_blue.png'
  },
  CIE: {
    icon: 'http://labs.google.com/ridefinder/images/mm_20_green.png'
  }
};

var cluster =[];
var min = .999;
var max = 1.001;

function load() {
  var map = new google.maps.Map(document.getElementById("map"), {
    center: new google.maps.LatLng(0,28.287167),
    zoom: 3,
    mapTypeId: 'roadmap'
  });

  var infoWindow = new google.maps.InfoWindow;
  var mcOptions = {gridSize: 20, maxZoom: 19};
  // this variable will collect the html which will eventually be placed in the side_bar
  var side_bar_html = "";

  // Change this depending on the name of your PHP file
  downloadUrl("phpsqlajax_genxml2.php", function(data) {
    var xml = data.responseXML;
    var markers = xml.documentElement.getElementsByTagName("marker");
    for (var i = 0; i < markers.length; i++) {
      var name = markers[i].getAttribute("assetno");
      var address = markers[i].getAttribute("address");
      var type = markers[i].getAttribute("type");
      var assetdesc = markers[i].getAttribute("assetdesc");
      var location = markers[i].getAttribute("location");
      var point = new google.maps.LatLng(
          parseFloat(markers[i].getAttribute("lat"))* (Math.random() * (max - min) + min) ,
          parseFloat(markers[i].getAttribute("lng"))* (Math.random() * (max - min) + min));
      var html = "
" + name + "
" + assetdesc + "
" + location + "
";
      var icon = customIcons[type] || {};

      setTimeout(function() {
        var marker = new google.maps.Marker({
            map: map,
            position: point,
            icon: icon.icon,
            animation: google.maps.Animation.DROP
        });
        return marker;

        setTimeout(function () {
            bindInfoWindow(marker, map, infoWindow, html);
            cluster.push(marker);
            // save the info we need to use later for the side_bar
            // add a line to the side_bar html
            side_bar_html += '' + name + ' - ' + assetdesc +'<\/a>
';
        }, 1500);
      }, i * 50);

    }
    var markerCluster = new MarkerClusterer(map, cluster, mcOptions);
    // put the assembled side_bar_html contents into the side_bar div
    document.getElementById("side_bar").innerHTML = side_bar_html;
  });
}

// This function picks up the click and opens the corresponding info window
function myclick(i) {
    google.maps.event.trigger(cluster[i], "click");
}

function bindInfoWindow(marker, map, infoWindow, html) {
  google.maps.event.addListener(marker, 'click', function() {
    infoWindow.setContent(html);
    infoWindow.open(map, marker);
  });
}

function downloadUrl(url, callback) {
  var request = window.ActiveXObject ?
      new ActiveXObject('Microsoft.XMLHTTP') :
      new XMLHttpRequest;

  request.onreadystatechange = function() {
    if (request.readyState == 4) {
      request.onreadystatechange = doNothing;
      callback(request, request.status);
    }
  };

  request.open('GET', url, true);
  request.send(null);
}

function doNothing() {}


texts:
I want to convert or a raster into a grid.
I am with QGIS. I have a raster of a parcel. I have a vector grid of the same area.
I want to add in each block of the grid a mean value of the raster.
Which tool or algorithm do i have to use ?  Maybe its cutting better than converting.


texts: Change from polygon to multipolygon in postgres

texts:
can somebody describe please how to do the transformation from polygon to multipolygon in the postgres DB, is not possible for me to save my edited data.

All data loaded just need to do some editing upon my polygons

Thanks in advance
regards


texts: Space-Time Principal Component Analysis with Missing Lat/Long Data

texts:
I am looking to run a space-time Principal Component Analysis on Shotspotter data from Brockton, MA: http://justicetechlab.org/data/. Shotspotter sensors record the timing, location, and number of gunshots. There are gaps in the data when no shots were detected.

I would like to include the gaps in my PCA. At the temporal scale of a day, I have filled in the gaps and assigned a value of zero rounds to days when there were no gunshots detected (e.g. row six of the table below).

I am confused about how to handle the missing lat/lon data for these zero shots-fired days. When there are no shots detected, there are no shots detected everywhere, so assigning a lat/lon seems strange. But I would like to include the zero shots detected data points.

Statistically speaking, would the best solution be to? :

 1. ignore the zero data completely
 2. use the mean lat/lon for the observed data
 3. randomly generate lat/lon data within the region of study and substitute these for missing lat/lon
 4. use another statistical method well suited to my problem

I have found some papers on how to impute missing measurements, but have not found much explaining how to handle missing lat/lon data:

http://pbil.univ-lyon1.fr/members/dray/files/articles/dray2015a.pdf
https://pdfs.semanticscholar.org/7bf7/e383847e7d29134f95894b463de731fd4c1d.pdf


texts: Fixing DXF import encoding problem using PyQGIS3

texts:
I have text with French accent in a DXF I import in QGIS3.16. The DXF is encoded in UTF-8 but as a result, the features attributes have encoding problems. The é charcacter are written Ã©.
Here's a portion of the code used:
# -*- coding: utf-8 -*-
layer = QgsVectorLayer(dxf_path + "|layername=entities|geometrytype=Point", "dxf_linestring", "ogr")
layer.setProviderEncoding(u'UTF-8')
QgsProject.instance().addMapLayer(layer)

When I convert the DXF from UTF-8 to ANSI in Notepad++, the problem is solved but I would prefer to convert in the script rather than converting the whole file.


texts: QGIS 3.6 picture in ui form (Qt5 Designer)

texts:
I'm facing trouble with display a picture in a ui form based on a field (path).
With the autogenerate form it works but when I try with a personalized form I cannot see the picture.
I tried with a Qlabel and QWebView in Qt Designer and in QGIS layer property I put the field as Attachment but without success.

I know that I can make a little python script to achieve this but I would like to store the ui in a db.

Does anyone have an idea ?


texts: Sharing model as toolbox in ArcGIS Desktop?

texts:
I created a model using ModelBuilder in ArcMap.

How can I share this model as a toolbox to run it in another computer?

The tool is expected to be shown under system toolboxes and users can run it as existing tools.


texts:
I need a sample map package that works with ArcGIS runtime 10.4.

I'm trying to develop a WPF/XAML application using the ArcGIS runtime SDK for .NET, but I've not been able to get very far because I don't have a map package to load into it.

At this point, I don't much care what map in particular, I simply need a map.

Our customer has sent us a .mxd with associated .gsd files, but when I try to create a map package out of it I get numerous errors: "00174 Label expression  uses VBScript, this is not support in the ArcGIS Runtime", etc.

I'd thought I'd be able to create a map in ArcGIS online to create a map package that would work with ArcGIS runtime, but when I open the map in ArcGIS desktop and try to create a map package, I get a different error: "00142 Layer's data source cannot be referenced in the ArcGIS Runtime".

Can anyone point me to a map package that will work with ArcGIS runtime 10.4? All I need is some street data, it doesn't really matter where.


texts: How to set up a querystring for a MapServer layer?

texts:
I am trying to query a MapServer layer via URL, but I am having trouble finding syntax examples. I want to use the indexquery mode and found this documentation http://gislars.github.io/mapserver_docs/cgi/controls.html as well as this http://mapserver.org/tutorial/example1-1.html. I am successful in returning the layer in map mode, but am not sure how to structure the indexquery string.

UPDATE: Upon further investigation it seems I need to use itemnquery mode. I have found some syntax that seems like it should work in this thread https://github.com/mapserver/mapserver/issues/4711, but I don't have access to the column names, and am not sure how to get them. Ultimately I'd like to query all records. Is there a general query for selecting all records, or is there a way to query a layer to return field names?


texts: Creating all necessary files from INTERLIS .ili file to produce and deliver data in QGIS?

texts:
I was provided with a INTERLIS 2.3 model as .ili file, and nothing else.

As far as I understand it, I should:

create a .imd meta file using ili2c.jar or the interlis Import/Export plugin by Pirmin Kalberer, i.e. the "Ili Model -> IlisMeta" tool
export an XML-file from the meta file (i.e. .xtf or .xml), encapsulating the data structure but no content yet using ili2pg.jar
import this structure into a postgreSQL / PostGIS database
create shapefiles (point, line, and polygon) using ogr2ogr.jar, which are (probably) already linked to the PostGIS DB
create the features (i.e., points, lines, polygons) and
use QGIS forms to fill the linked database with attributes

I already fail at the first point. Java delivers a .NullPointerException, and the plugin doesn't create any output without a warning or an error.

How do I produce the necessary files properly?

I'm completely new to INTERLIS, and might have misunderstood some concepts, but would be happy if someone could point me to the right ressources to get to productive GIS layers from an .ili file. I'm running QGIS 2.14.8-Essen (LTR), and the INTERLIS model in question is not (yet) published.


texts:
In QGIS 3.6.2 (on Windows 10), one can open the properties of a layer, click on Metadata, and enter information about the layer. At the bottom of the properties window is the option to save metadata as a file (a QMD file) or "Save as Default".

Where are default metadata stored for PostGIS layers in @qgis 3.x? It's not in the project file because I didn't save a project file but when reloading layers from PostGIS, the metadata changes are preserved.


texts: ArcObjects - cheapest way to check if license has been checked out

texts:
So I noticed a big performance hit when I add the following line to my code:

ESRI.ArcGIS.RuntimeManager.BindLicense(ESRI.ArcGIS.ProductCode.EngineOrDesktop);

I have two separate classes that call the BindLicense method.  After the first time it is called, it isn't necessary to call again, but neither class can assume that has happened.  What I want to know is if there is a cheaper way of checking if a license exists before attempting to bind?

I feel like I've seen code to do this in Python, but don't know if there is the equivalent in ArcObjects...


texts: How to convert a "yyyy/mm" field to "yyyy/mm/dd"?

texts:
I am trying to edit an existing date field of the form "yyyy/mm" to a format that the TimeManager Plugin for QGIS can handle - I guess "yyyy/mm/dd" will do it.

I've got as far as opening the attribute field table, adding a new table and getting the field calculator up - but I can't see how to simply add a "/01" to the existing info in the field without getting a "invalid expression" error.


texts: Broken MSSQL provider in recent QGIS versions - has anyone experienced this?

texts:
I recently moved from QGIS 2.2 to 2.4.0 and in the process the MSSQL spatial data provider appears to have stopped working properly (behaviour is also the same with 2.5.0).

QGIS will still connect to an MSSQL table and will display the geometry, however it is not possible to filter it with a query, the SRID is not being read from the geometry_columns table and more importantly it is not possible to commit any edits back to the database.

I am running SQL Server 2012 - the geometry_columns table is present and populated correctly and all tables have integer primary key columns, etc.

Has anyone else encountered this problem? I am a bit perplexed as I have never experienced any problems prior to updating to 2.4.0.


texts: How can I fix the ArcToolbox error TypeError: cannot concatenate 'str' and 'NoneType' objects?

texts:
I am trying to run the ExportHTMLImageMaps tool in ArcGIS 10.3.1 to create HTML image maps to be used as static web geographic navigation maps, but I keep getting the Error TypeError: cannot concatenate 'str' and 'NoneType' objects

I checked each field being used in the tool and each is type String, so I'm not sure where to go from here.  Below is a screen cap of my error log.  I am a novice at Python and troubleshooting code.


texts: Using initQgis() on headless installation of QGIS 3?

texts:
When running PyQGIS testing 2.99 the following coredumps:

from qgis.core import *
QgsApplication.setPrefixPath("/usr", False)
qgs = QgsApplication([], False)
qgs.initQgis()

How to start migrating PyQGIS scripts to the upcoming QGIS 3?


texts: Slow execution of request for searching nearest street to GPS point

texts:
I use OSM data. I need to find nearest street for GPS point. I wrote next SQL query:

SELECT
     osm_id,
     name,
     ST_DISTANCE(geometry, ST_SetSRID(ST_POINT(37.72308, 55.47957), 4326), true) as dist
FROM
     roads
WHERE
     ST_DWithin(geometry, ST_SetSRID(ST_POINT(37.72308, 55.47957), 4326), 1) ORDER BY dist LIMIT 1;

But it's execution is very slow. It's take 7.5 seconds for not big territory. Table structure is next:

newdb=# \d roads
                                   Table "public.roads"
  Column  |           Type            |                     Modifiers
----------+---------------------------+----------------------------------------------------
 id       | integer                   | not null default nextval('roads_id_seq'::regclass)
 osm_id   | bigint                    |
 type     | character varying         |
 name     | character varying         |
 tunnel   | smallint                  |
 bridge   | smallint                  |
 oneway   | smallint                  |
 ref      | character varying         |
 z_order  | integer                   |
 access   | character varying         |
 service  | character varying         |
 class    | character varying         |
 geometry | geometry(LineString,4326) |
Indexes:
    "roads_pkey" PRIMARY KEY, btree (id)
    "roads_geography" gist (geometry)
    "roads_geom" gist (geometry)

is there any way to increase speed of the query?

PostGIS 2.2
PostgreSQL 9.6


texts: Reclassifying Raster (Heatmap) from Value to Percent using QGIS?

texts:
I have a point layer with several points "belonging" to different years by attribute. For all points of each year I created a heatmap. Now I would like to define the core area of each heatmap.

Unfortunately, the values of the different heatmaps are pretty different, therefore its not possible to choose one value for all. But because there are over 100 different years, so over 100 different heatmaps, I do not really want to set the value for each core area manually.

I thought of reclassify the heatmaps from value to percentage (min value = 0%; max value = 100% or min=0 max=100) to generate the core area based on one percent-value for all heatmaps with batch-processing or graphical modeller.

I tried different recalculating/reclassifing-tools like the raster calculator, r.reclass (grass), reclassify values (saga) and reclassify by table (qgis) but I'm not able to get it to work. But honestly I think I'm just too stupid to find out the respective expression for this problem.

Can anybody help me with the correct expression or way to reclassify the rasters/heatmaps from min-value=0(%) to max-value=100(%)?

I prefer working with the new QGIS Version 3.6 but if a solution requires another version of QGIS it would not be a problem.


texts:
I met a problem when I was using this equation f=(NDVI-NDVI MIN)/(NDVI MAX-NDVI MIN) to get vegetation cover.
The problem is “NDVI-NDVI MIN” part. NDVI value range is (-0.05115,0.605945)When I used Raster Calculate to run it, the result should be (0,0.657098), but it became (-1,49012e-08,0.605945).
It looks like there is some problem when NDVI value become 0
Any ideas about which part I did wrong?


texts: Difference between node and vertex in GIS?

texts:
In graph theory the term of node and vertex are equivalent, but sometimes in the GIS world the nodes and vertex are mentioned in the same context.

I was wondering what the difference is?

Node: can be a start and an end node from a link (edge)


texts:
I am wondering if there is a way to convert a graphic produced from symbology in QGIS to a vector layer (ie: shapefile).  I have a point layer with bearing direction in the attribute table.  Following the answer from How to display an arrow pointing to a specific bearing in QGIS, I have been able to make a line graphic depicting the look direction from a photo collection.   I would like the line graphics to be converted into a vector layer, so I can perform some spatial analysis.

Any thoughts?

Thanks....


texts: Multiply two image collections with different sizes based on dates

texts:
I have two image collections. They are both filtered by date beginning on Jan 1, 2019, and ending on today's date. The size of the first (kc) is 35 images and the second (ETo) is 225. The first has an image every 16 days, whereas the second collection has one daily.

I am trying to multiply in a loop the last kc image by groups of 5 ETo images by looking at their date. If the current ETo image in loop has a date more recent than the last kc image, move to the next kc image and wait for a new increment. This is the logic I thought of:

newCollection = []
i = 0
currentDate = kc[0].date
foreach ETo as etoImg:
    if (etoImg.date >= kc[i].date):
        i++
        currentDate = etoImg[i].date
    newCollection.push(etoImg * kc[i])

Having said that, this is the code I made for calculating that:

var kclist = kcAccum.sort('system:time_start').toList(kcAccum.size());
var kcindex = 0;
var ETaColl = ETo.sort('system:time_start').map(function(etoimage) {
  var currkcimage = ee.Image(kclist.get(kcindex));
  if(currkcimage.get("system:time_start") < etoimage.get("system:time_start")) {
    kcindex++;
  }
  return (etoimage.multiply(currkcimage))
    .copyProperties(etoimage, ["system:time_start"]);
});

var ETaChart = ui.Chart.image.series(ETaColl, geom);
print(ETaChart);

The problem is that once the graph is printed, it shows very strange values, very different from what was expected.

The first value shown in the ETo graph is 1.953. The last in kc is 0.228. Therefore, the first value in the resulting graph should be 0.445. However, it has a value of 0.374. Similarly, the last value shown in the ETo graph is 6. The last in kc is 0.103. Therefore, the last value in the resulting graph should be 0.618. However, it has a value of 1.149.

That is how I concluded the calculations were wrong.
Does anyone have any suggestions as to what I may be doing wrong here?

I have heard it is possible interpolate image collections to generate images for the dates where the kc collection has no images. However, I don't really know how I could do that. I am open to suggestions.

This is the ETo graph:

This is the kc graph:

resulting:


texts: Dissolve polygons with road network gaps

texts:
I am using QGIS and have a shapefile of the blocks in a large city, where each polygon is a particular streetblock. The city is divided into various sectors, and each block has information on its sector. I can easily classify the blocks by this attribute to visually inspect the shape of these sectors:

The issue is that the gaps between the blocks is the road network. I would like to dissolve all blocks in the same sector into a contiguous polygon, but am unsure how to do this.

Does anyone know an easy way to "fill in" the gaps between the city blocks in order to make the dissolve work properly?


texts: Overstory-understory classification from discrete return aerial LiDAR

texts: I have an opportunity to collaborate with bat biologists who are interested in predicting probability of roost occurence in seasonally flooded ponderosa woodland. We suspect that both properties of canopy and understory are relevant, and I hoped to extract metrics from each, eg, rumpleness and density, from a discrete return aerial LiDAR point cloud. My impression of the system is that there are single overstory trees with a predominantly herbaceous understory, and I'm wondering if there are useful tools for binarizing the cloud into overstory and understory in such a system. I've found recent methods that look for dips in histograms of Z values in a moving (XY) window, and make a cut at the dip. I'm wondering if there is a good candidate algorithm (currently published and implemented) that could perform this task?


texts: Assigning RGB values from Geotiff image to LiDAR data, using R

texts:
I have given a Geotiff image and its corresponding Lidar data (x,y,z) in UTM coordinates. I need to merge the Lidar data with the RGB values from the image.

That means, at the end, I need to plot (3D) each point of the LiDAR cloud color coded with its corresponding RGB value from the Geotiff image.

I converted the Lidar data into a shapefile using QGIS. What should I do next?

In R, I tried the plot3D function, but, it did not work. I'm attaching the text doc, shapefile, and tif image

Edit:

I've done the following program as shown below:

require(raster)
require(maptools)  # to take shape files
#require(car) # for scatter3D
require(plot3Drgl)

##setwd("C:\\Users\\Bibin Wilson\\Documents\\R")
##source('Lidar.r')

data = read.csv("C:\\Users\\Bibin Wilson\\Desktop\\Lidar\\lidardata.csv")
#nr = nrow(data)
nc = ncol(data)

nr = 500

require(rgdal)
X = readGDAL("C:\\Users\\Bibin Wilson\\Desktop\\Lidar\\image.tif")

topx = 4.968622208855732e+05;
topy = 5.419739403811632e+06;

final = matrix(nrow = nr, ncol = nc+2)

for(i in 1:nr) {
 x = data[i,1]
 y = data[i,2]
 rr = round((topy-y)/0.0833)
 cc = abs(round((x-topx)/0.0833))
 if(rr == 0) {
  rr = 1
 }
 if(cc == 0) {
  cc = 1
 }
 final[i,1] = x
 final[i,2] = y
 final[i,3] = data[i,3]
 final[i,4] = rr
 final[i,5] = cc
}

for(i in 1:nr) {
 x = final[i,1]
 y = final[i,2]
 z = final[i,3]
 rr = final[i,4]
 cc = final[i,5]
 if(rr <= 5086 && cc<=3265) {
  r = X[rr,cc,1]/255
  g = X[rr,cc,2]/255
  b = X[rr,cc,3]/255
  c = cbind(r,g,b)
  scatter3D(x,y,z,2,c)
 }
}

But while trying to plot the graph, it shows the following error:

  Error in [.data.frame(x@data, i, j, ..., drop = FALSE) :
    unused argument (1)

Edit:

I got the 3D model without the RGB as shown below:


texts:
I am trying to shuffle data in a file with four columns as follows:

A  1  2  3

B  0  5  6

C  2  4  7

D  7  8  9

E  10 4  8

I only want to shuffle the first column. So, am expecting something like this

E  1  2  3

D  0  5  6

C  2  4  7

A  7  8  9

B  10 4  8

Does anyone have any idea how to do this? I have tried something like this

with open("file.txt", "rb") as source:
    lines = [line for line in source]

import random
random_choice = random.sample(lines, 5)

with open("newfile.txt", "wb") as sink:
    sink.write("\n".join(random_choice))

with open('newfile.txt', 'r') as fd:
    for l in fd:
        if l != '\n':
            print(l.strip())

But observed that all the other columns also changes.


texts:
I am trying to match a large set (43M) of small circular polygons to a smaller set of very large polygons.

To make this process as fast as possible I've created points from the centers of the small polygons and have successfully matched them to the the larger polygons substantially faster than before using the following code:

Large_poly_layer > The very large polygons
pts > points data source

for large_poly in Large_poly_layer:
    polygeom = large_poly.GetGeometryRef()
    pts.SetSpatialFilter(polygeom)
    for pt in pts:
        fid = pt.GetFID()
        // this then writes the FID to a numpy array
pts.SetSpatialFilter(None)

Now what I am trying to do is obtain any points that didn't match, so that I can buffer them into polygons and see if they'll match at that points.

I think the fastest way to do this would be to do the inverse of the spatial filter but I cannot find a decent way to do it.

I have also attempted the following:

-copying the whole large point dataset to a memory layer in ogr then deleting points that have matched, but the lack of spatial indexing in the memory format makes this process incredibly slow.

-I also cannot write back to the original source data as it is contained in a .gdb file


texts:
I am trying to create a cloud free image to use in the unsupervised classification, however my final product ("composite collection without cloud") stays with stranger stains and it seems to influence on the classification.
Does anyone know the reason for the stains, and how to soften them?

Link to my code:

https://code.earthengine.google.com/c18d6de5ada3e21a5b1aca9786efb103


texts: Joining/Extracting Attribute Data from two SpatialLinesDataFrames with no common Attribute

texts:
I have two rather large datasets, but I will only show a "single" subsample here.
I have one SpatialLinesDataFrame (RFI) shown here:

and another (CLASS):

These two line features are the same section of road. The CLASS though is broken down into 15 segments along that road and the RFI contains other important information, like the road name.
My initial thought was that it would be simple to just perform an intersect on the two datasets, so I first changed the column names in the RFI so they would not conflict with the matching ones in CLASS.
#CHANGE SOME HEADER NAMES SO THEY DO NOT MATCH MAINTENANCE CLASS DATA
names(mot_rfi)[13] <- "RFI_HWY_NUMBER"
names(mot_rfi)[14] <- "RFI_OBJECTID"
names(mot_rfi)[15] <- "RFI_FEAT_LEN"

I then learned that intersect does not work, which lead me to gIntersect, which produced and error
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : spgeom1 and spgeom2 have different proj4 strings
As my proj4 strings appear to match I fumbled around for a quite, then decided that maybe a SpatialPolygonDataFrame would work better as many examples use a SpatialPointsDataFrames and SpatialPolygonDataFrames. I fumbled around for a while and so far the closest I have come to my intended result is:
sf_CLASS <- st_as_sf(CLASS)
sf_CLASS_polygons <- st_polygonize(sf_CLASS)
sf_RFI <- st_as_sf(RFI)
sf_RFI_polygons <- st_polygonize(sf_RFI)
G <- st_join(sf_CLASS_polygons, sf_RFI_polygons, left = TRUE)

Where my DataFrame has the original full CLASS dataset, plus the headers for the RFI dataset with no values!
I also just noticed I failed to change one matching column header, but R did it for my by adding a .x, and .y.  Ill have to go address that in another names line.
So why is it that I can not find a way to accomplish this task? It seems so simple in my mind. Extract data from RFI into Class based on geographic location, not similar attribute. What is the simple solution I have failed to uncover? Is a Line over Line not a common task?


texts: QGIS and SQL Server Spatial Tables Performance

texts:
I'm trying to get a SQL Server 2008 R2 spatial database up and running to streamline the production of maps using QGIS.

Currently I just have a handful of spatial tables (a postcode level point table and a few polygon tables for various boundary levels).

My initial plan was to use these to create views (primarily point in polygon counts based on UK postcodes) which I could then pull into QGIS to make the data processing side of map production a lot faster and less fiddly.

Unfortunately my plan has been scuppered somewhat by the amount of time it takes just to connect to my SQL Server using the MSSQL layer tool. QGIS seems to process each table/view prior to allowing me to select any of them. Not only is this extremely slow, it also means that each time I add a new view or table to my database, the whole process is going to get slower and slower until it becomes unusable.

Does anyone have any pointers on how I might be able to speed the QGIS-SQL Server connection process up?

I am beginning to think that saving views is a bad idea and that I should maybe use 'select ... into...' scripts to generate tables and minimise the processing required. I'd be interested in peoples' views on this.


texts: Measure the area of a polygon resulting from the intersection of two other polygons?

texts: How can I measure the area of a polygon resulting from the intersection of two other polygons? If I use the measure tool or open the attribute table for the resulting polygon, it shows me the areas of the original two polygons...


texts: Combining parcels that are adjacent and add up to more than 10 acres

texts:
I am using ArcMap 10.7.
My task is to find a site to be developed for parking. Sites are supposed to be greater than 10 acres. But we can also look at adjacent sites that add up to 10 acres. I read something similar in this link:
Merging adjacent polygons to one which sum of attribute values of adjacent ones is between 1000 and 2000?
However, I have not been sure how to apply it. I have tried another rudimentary way which is to select parcels between 0 and 1 acre, then select parcels next to them, then select from them those which are 9+ acres. Then repeat the process from 1 to 2 acres, 2 to 3 acres, etc. However this is time consuming and this leaves a little room for error if say I get a parcel of .05 acres next to one that is 9 acres and they don't add up to more than 10.
I would love to just be able to do this all from the parcels layer. But if there's a way to create the "0 to 1 acre" layer and then select from the parcels layer something like "0 to acre parcels" acres + "parcels" acres >= 10"


texts: Testing for polygon orientation using GeoTools

texts:
I want to merge multiple geometries with Geotools (starting off with one geometry) but the result is a polygon that has a clockwise orientation. This gives an invalid geometry.
This is my code for the merge:
    GeometryCollection geometryCollection =
    (GeometryCollection) factory.createGeometryCollection(list.toArray(new Geometry[list.size()]));
    Geometry all =  geometryCollection.union();

Why would something give an invalid geometry and what can I use to test for orientation?


texts: Appending country names to route feature class using ArcGIS Desktop?

texts:
I have a polyline feature class made up of multiple segments. This feature class overlays a polygon feature class representing the countries through which this polyline is going. For each segment, I would like to add an attribute containing all the countries through which the segment crosses, delimited by a comma for example.

I have tried a Spatial Join, using the polyline as the target and the polygon as the join feature. Unfortunately this has resulted in null values for the country attribute.

Any ideas?


texts:
Is there a way or a tool to mimic "find and replace" functionality in a Geodatabase (File, Persoanl, Enterprise). For example I need to remove a string from all FeaturClass names.

for example in the following cases, I need to remove "Pattern_":

Pattern_Road1
Pattern_Parcel
Pattern_Buildings


texts: Is there an equivalent of arcgis' eliminate feature for QGIS?

texts:
I already know of grass' v.clean, but I'm looking for something that doesn't involve third party software.

Python scripts or addons are good candidates.

What is arcgis eliminate : http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Eliminate%20(coverage)


texts: Commandline tool to change RGB value of GeoTIFF

texts:
I would like to change the rgb-value of raster-tiffs from 0,0,0 to 1,1,1.
Is there a commandline-tool that can change the rgb-values of GeoTIFF files?

At the moment I use irfanview to change the color. As long as I keep the TFW-file of the GeoTIFF but do change the color for many raster files is not very conveinient.


texts: Putting multiple condition with more raster to Con tool in ArcGIS Raster Calculator?

texts:
I have three rasters and need to make a condition including all of them. The condition should look like Ras1>(Ras2*85 + (Ras3/4)), Ras1 if true, 0 if false.

Can anybody here help me to do this?


texts: Error using calc() function: cannot use this function using R?

texts:
I need to run calculations on large multi-band rasters and export a RasterBrick, and am trying to do so using the calc() function in the raster package, for the purpose of memory efficiency. The function runs fine on its own, but when I try to include it in calc(), I keep getting this error:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function

How can I make this work?

Simplified code:

fn = system.file("external/test.grd", package="raster")
s = stack(fn, fn, fn, fn)

out = calc(s, fun = function(x){
  for (i in 1:nlayers(x)){
    x[[i]] = x[[i]] - cellStats(x[[i]], "min")
    x[[i]] = x[[i]]* 5
  }
  list = unstack(x)
  out = brick(list)
  return(out)
}
)

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function


texts: Boundaries of Northern Ireland Local Government Districts (LGD2014)?

texts:
I'm looking for boundaries for Northern Ireland's Local Government Districts (LGD2014).
I only find a list of units, but no the boundaries. There are 11 LGDs.

Does anybody know where to find the boundaries at a decent scale for analysis?


texts: Given (lon,lat), how can I quickly retrieve the containing administrative region?

texts:
Given a coordinate (lon,lat), I'd like to retrieve the containing administrative region (state/province level) for an offline mobile application. The method employed needs to be fast/efficient. NaturalEarthData is one source I found that has the required shape/attribute data (http://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-1-states-provinces/).

I'm looking for suggestions on how to convert (lon,lat) into its corresponding administrative region. The obvious way is running point in polygon and maybe preempting it with a bounding box check. Are there better ways?

edit:

Spacedman:
When I meant mobile I meant device side; right now I'm targeting embedded Linux, Android and QNX (BlackBerry). I think GEOS might be overkill -- I don't know if its worth adding all of GEOS for spatial indexing. In either case, right now I'm trying out this method and R.K.'s raster approach as well.

R.K.:
This seems like a great idea, especially since the data is static... but I'm having a bunch of trouble trying to generate an image to do the lookup. I'm a novice with foss gis tools but I haven't been able to generate what I need with Mapnik or QGIS so far.

The NaturalEarthData shapefile separates each admin region into separate geometries with unique ids. I added an attribute column to this shapefile that converts the uid directly into a color, thinking it'd be possible to render a map where I could set the color of each polygon based on the attribute... that assumption didn't pan out too well. I'll keep at it though, and update this question with a properly selected answer soon.

edit2:

Alright, so with the preliminary testing I did, I think that the raster approach worked out a bit better for me. The size of the data ended up being much smaller compared to the vector data and the look up times were pretty quick (most of the time is spent loading a tile into memory). The overall complexity of the raster implementation is pretty low as well. The time taken for point in polygon tests for vector data varies a lot based on polygon shape, and you need to watch out for multipolygon definitions / erroneous geometry as well.

I created a really large map (36000x18000) so I got an acceptable level of accuracy (100px/degree latitude) with respect to the consideration of geopolitical borders. The raster data + meta data (admin region name) only took up 2.5mb.


texts: Getting the feature information under the mouse (GeoTools)

texts:
I'm interested in programmatically finding out information of the features (of a map) under the mouse. I'm using GeoTools for drawing the map using a shapefile. I've first used the following code line by overriding onReportedUpdated().

InfoTool tool = new InfoTool() {
        @Override
        public void onReporterUpdated() {
            String text = getTextReporterConnection().getText();
            System.out.println(text);
        }
    };

But, here I get a full string with all information. Yet, I would like to individually find those feature information.

I've found this question, in which it is said that there's no direct API call for this. But since this is a quite old question, I would like to know whether there's any direct API methods to achieve this.


texts: How to quickly extend PostGIS's R-Tree implementation?

texts:
I am working on a problem where along with data points location, I also have information about the connectivity between points.(Think a graph of points, where an edge denotes some relationship)

My queries have to be optimized not just to find say 'Nearest Neighbors(NN)' or 'Reverse Nearest Neighbors(RNN)' but to find the points that satisfy NN/RNN along with the graph constraints.

It seems like I might need to modify the R-Tree index to store this additional information and perhaps even write a couple of my own routines to add some extra functionality.

Is there any way I can abstract the implementation of R-Tree by PostGIS and quickly extend it? Do I need to delve into the internals before I can get anything accomplished?


texts: How can I create a QGIS plugin that needs an external Python library?

texts:
I’m trying to enhance QGIS's Profiletool plugin by adding DXF export functionnality. For this, I need to to import dxfwrite library.
As it's installed on my computer, i added the following lines to plugin code:

import sys
sys.path.append("C:\Python27\Lib\site-packages")

With those short lines of code, I'm able to load dxfwrite. But it's hardcoded, so if the user does not have Python 2.7 but 3.4 (for example), library loading will probably fail because plugin will try to load dxfwrite from a non-existing repository. The same problem may happen if python path is not C:\Python27\....

I tried to put dxfwrite sources into plugin directory, but 1) with no results, 2) if results, it would link plugin to a specific version of dxfwrite and 3) if results, it would increase a lot plugin size.

So, how can I properly load dxfwrite without facing those problems ?


texts: Speed- and COG-aware trajectory interpolation/resampling using Python

texts:
I have a set of ship trajectory points that were irregularly sampled and I would like to resample them to a fixed frequency (e.g. every minute).

Rather than just linearly interpolating between latitude and longitude, I would like to do this in a way that smooths according to course (COG) and speed (SOG).

Here is an example of the data:

Time              Lat    Lon    COG  SOG
2019-05-08 00:16  36.25  14.82  211  8.0
2019-05-08 00:28  36.23  14.80  221  7.7
2019-05-08 00:41  36.21  14.77  222  8.3

And here is a toy diagram explaining what I am trying to do. Instead of the first, I would like the second:

Is there a simple way to achieve this using Python libraries? I am hoping there is a built-in method out there that can do this in a few lines :)


texts: I downloaded an SRTM DEM GeoTIFF file from Earth Explorer. It looks fine when previewing it on Earth Explorer (see picture) but when added in ArcGIS Pro the DEM looks too bright in some of its tiles.


texts: Openlayers 2 projection different depending on how I ask for it

texts:
I have a map and am trying to see what projection is set, but, depending on how I ask, I get two different values, this is what I see in the console:

> this.olMap.projection
< "EPSG:4326"

> this.olMap.getProjection()
< "EPSG:900913"

Is this expected? Which one is true?


texts: Is it possible to rename a field with the current date

texts:
I'm trying to figure out if it is possible to rename a field with the current date.  What I am trying to accomplish is a table to table conversion, dropping a dozen fields, leaving only two fields, Name and Status.  I want the Status field to be renamed to Status_yyyy_mm_dd.  The script will then run a join field to another table, thus creating a history of the status field by day.  I am running Arc 10.1.  I have everything working, except the naming bit.  Any help would be greatly appreciated.

import arcpy
import time

uiq = "uiq_data_set"
gdb = "c:\\Temp\\default1.gdb"
newFieldName = "Status_" + time.strftime('%Y_%m_%d')

arcpy.TableToTable_conversion(uiq, gdb, "uiq_historic", "", "Name \"Name\" true true false 255 Text 0 0 ,First,#,c:\\Temp\\default1.gdb\\uiq_data_set,Name,-1,-1;Status\"Status\" true true false 255 Text 0 0 ,First,#,c:\\Temp\\default1.gdb\\uiq_data_set,Status,-1,-1", "")


texts: Error while displaying qgis project via qgis2leaf plugin?

texts:
couldn't see the map on the browser while using QGIS2LEAF Plugin in QGIS 2.10 Pisa. It gives error on web page saying:

  error at line 8: L is undefined.

Please let me know if anyone knows the solution.


texts: Filtering dataset using SQL or Python in ArcGIS Desktop?

texts:
Using ArcGIS Desktop 10.0 SP4.

I have a large dataset of sites that needs to be filtered down to a usable size.  There is an attribute for the site_ID and another one for the site_type.  There are duplicate site_ID’s in the dataset and I only need one record for each site.  The preference of what gets deleted depends on the type.  If there is a record with a type “A” then all other records for that site can be deleted.  If a site does not have a type “A” then it doesn’t matter which records get deleted.

In a nutshell I want to turn this…

Site_ID-----Site_Type
  123--------------A
  123--------------B
  123--------------C
  555--------------B
  555--------------B
  555--------------C

Into this…

Site_ID-----Site_Type
  123--------------A
  555--------------doesn’t matter

I’ve been playing around with SQL statements without any luck.  Any suggestions?


texts: openlayers 3 (or oracle) convert 2d to 3d

texts:
We have an existing oracle gis which is being reworked. We're not really touching the oracle stored procedures or tables.
In the table, the sdo_geometry appears to be a 3D line, with 0 as the z value. I'm not entirely sure why it was written this way, I wasn't there when it was written.  However, when this data gets brought into an openlayers vector layer and I observer the feature geometry, there is no Z value, so somehow it's being converted from 3D to 2D.
In my new code I'm actually sending the data to oracle as a WKT, and then converting it to an sdo_geometry in an oracle helper function which passes the data to the other procedure which takes an sdo_geometry (I couldn't find a way to convert to sdo_geometry format in c# or javascript). So the linestring is being sent as 2D which the database expects 3D so therein lies my problem. I think I can get around this on a few ways, but I don't know how to do any of them. I could...

A. Force my vector layer in openlayers to be 3D. with z value of 0.

B. use some (unknown to me) openlayers function to return a 3D version of the geometry (or WKT) given a 2D, with Z set to 0.

C. use Oracle (in my helper function) to make the sdo_geometry 3D (a different constuctor maybe?) currently using sdo_geometry(wkt_in,srid) srid is 8307 if that matters.

D. Write my own function in javascript or C# which will just parse the WKT and add a 0 to each coordinate set. I don't have a problem doing this, but would rather use existing functions if they're available.

E. something else I haven't thought of :)


texts: Change SDE Version in ArcGIS JavaScipt API depending on connected user

texts:
I am programming for an web editing application with ArcGIS API for JavaScript. I've published a version of ArcSDE with Feature Acces and I need to know how can I change in the web application from this version to another through the user that is connected to the application using the security manager of ArcGIS Server.

I've seen that exist a parameter gdbVersion and SetGDBVersion of featurelayer, but I don't know how to use. I need to see some example of this because the information in the documentation is not enough.

Can anyone help me? Any sample, please?
Thanks


texts: Batch select by attribute and batch export output

texts: I have ArcMap and a shapefile attribute table with lots of columns. Each column contains values between -1 and 1. For each column, I would like to select the values greater than 0 and then export a new shapefile from that selection. I have about 100 columns so don't want to do it one by one. Do I have to use ModelBuilder or python to do this?


texts: What is typical v.surf.xxx execution time?

texts:
What is the typical v.surf.xxx exicution time?

I am trying to run v.surf.bspline on a set of ~1600 points.  It has bee running for about an hour now and is stuck at spot indicated in the image.

I am wondering about the time because I am not sure which interpolation method I will use and want to see the results fom all of the diferent methods.


texts: Leaflet kmz stops loading other kmz files if gets XHR 404 while fetching remote file in any one of kmz's

texts:
Referring to https://github.com/Raruto/leaflet-kmz
Here is my leafllet kmz code:

Referring to the code above:
kmz.load('https://raruto.github.io/leaflet-kmz/examples/blah.kmz');
ends in 404 error and stops subsequent loading of kmz.
kmz.load('https://raruto.github.io/leaflet-kmz/examples/globe.kmz') doesn't get loaded because of 404 error while loading blah.kmz file. I don't want the code to break in case of 404 error in any of the kmz files. I want to ignore 404 error and proceed to loading of next kmz.
How that can be implemented here?


texts: Tkinter crashing ArcMap from Python Add-In

texts:
ArcMap crashes every time I try to display a form using tkinter from my python add-in button.

Are there any updated workarounds for this?

I saw this from before, but couldn't get it to work:
https://gis.stackexchange.com/a/212848/108211

Has anyone figured out a way to make this work, or are there any other python GUI's that will work with python addins?


texts: Toggle Layers in OpenLayers does not work?

texts:
I am trying to make a map in which I can toggle layers on and off in the navigation bar. A simple version of my code is in the following JSFiddle: http://jsfiddle.net/TimLucas/3cqsmc2L/

If you click on the button in the nav-bar, the layer should disappear, but for some reason it does not. I am using this type of nav bar since I would like to have sub sub sub menus (this site is not to be used on mobiles).

For some reason I cannot figure out why the toggling does not work, I've been looking for 2 hours now but I cannot find an awnser. The toggeling of the layers does work in the console (if I run the file locallly) so this is probably related to the eventhandler on the navbar.

Can anyone see what the problem is?


texts: Why does quantities graduated colors shows different data?

texts:
I am currently working on a heatmap and trying different ways of visualizing my data. I used unique values in combination with labels which worked fine.

However, it may be to confusing for a research paper and and I decided to use quantities instead. Unfortunately, using the same value I get very strange numbers that are not present in my dataset.

I haven't used ArcMap in while but this does not seem right to me. Does anyone know what to do?


texts: Manage attribute fields ArcGIS Pro 2.0.1

texts:
I'd like to edit, create new and delete attribute fields using ArcGIS Pro.
I open attribute table using short cut Ctrl + M.
When I click Field View

So I try to add new, or any manipulation. But the most important question. How to save edits?
ESRI tell me to click save. But where's "save" button?


texts:
Currently I am working on a project in Kosovo. After gaining autonomy the country "invented" their own CRS, which is driving me crazy sometimes. It doesn't even have an EPSG-Code, it goes by "Kosova Ref 01" and the SR-ORG 7392, 9079 and 9080.

I managed usual reprojecting, but today I encountered a significant difference in altitude. While e.g. Google Earth and some surveyor-data state, that the area is I am working on is between 1450 and 1600 m above sealevel, data using the Kosova Ref 01 says it is about 40 to 45 m higher. Is this possible, or did someone (maybe me) botch something?

tl;dr: Is an altitude difference of about 40 m between two CRS feasible?


texts: Equivalent of CORINE land cover data sets for rest of world

texts: I am looking for the equivalent of CORINE land cover data sets (100 meters) for the rest of the world; not merely urban areas but any area in general. I see sousa's question unanswered so I thought I will ask again.


texts:
When I run the field calculator to update area, the resulting areas are too large. This happens whether I use QGIS 2.18 or QGIS 3.4.

The default CRS is set to: ESPG:27700, OSGB 1936 / British National Grid.
However, somehow, the CRS has changed to ESPG:4325, WGS 84. When I change it back, the polygons do not appear on screen. If I recalculate $area, the amounts are '0.0000'.

Area default set as hectare, layer is a vector layer and all polygons are single part.

Example: using measuring tool gives area of 1 polygon at 7.441, but field calculator $area = 12.1819.

Any suggestions?


texts:
i have wms layer with coming geoserver and it's bounds is

 bounds = new OpenLayers.Bounds(417138.6205791628, 4243634.852681775, 486211.1950145211, 4331542.555877601).transform(WGS84, WGS84_google_mercator);

but when i use google baselayer, in the opening map it shows me full extends how can i transform my bound with harmonic google. ?


texts: GeoServer function PointN doesn't work in YSLD

texts:
I did some tests and want to confirm if function PointN doesn't work in YSLD. For example below YSLD doesn't display any label:

name: style
feature-styles:
- name: style
  rules:
  - name: style
    symbolizers:
    - text:
        uom: metre
        geometry: ${pointN(geometry,1)}
        label: ${geometryType(geometry)}
        placement: point
        fill-color: '#000088'
        font-family: Tahoma
        font-size: 0.5

Change the geometry and label part to display geometry type also won't display any label:

geometry: ${geometry}
label: ${toWKT(pointN(geometry,1))}

But changing label to ${getX(pointN(geometry,1))} will display always 0.0 regardless of number of point we enter:

geometry: ${geometry}
label: ${getX(pointN(geometry,1))}

I was looking already for YSLD binaries but I've found only GeoServer binaries.


texts: How to generate a layer of points from a database of lat, long using qGIS?

texts:
I am a programmer but I am new in QGIS. I have a database of some hundreds of points with lat, long information. It is an excel spreadsheet.

How can I create a layer of points of each lat, long from those data?


texts:
I'm trying to create an Aspect Map similar to this one:  https://kingsgeocomputation.files.wordpress.com/2016/03/10sikkim-village-aspect.jpg?w=980

I have two files to work with from the start; the primary is a contour .dwg that I know is accurate while the secondary is a DEM that was created in Civil 3D using those contour lines (Civil made a huge file out of it so I'm not a big fan of using it if I can start from scratch using the contours in QGIS).

In my first attempt I used the DEM and the Aspect tool within DEM (Terrain Models) to create a map.  It worked, but for some reason it looks like it's converting the DEM back to contours and coloring the contours rather than the entire surface (as if it interprets the slope between the contours to be flat).  This is a screen capture of what that's looking like:  http://i535.photobucket.com/albums/ee355/MR_HT/Aspect%20Map%20Attempt%201_zpsrikymkrj.jpg

It works in the steep areas, but the flatter areas don't really cut it.  This is a screen cap of the DEM I used to make this:  http://i535.photobucket.com/albums/ee355/MR_HT/Aspect%20Map%20DEM_zpsuushwzip.jpg

I presume that the problem is that the DEM was generated as flat tiers between the topo lines.  So it looks pretty good as a greyscale but when QGIS tries to interpret the slopes it's reading each grey as a flat plane instead of smoothing the colors across the surface.  Does this make sense?  If this is a problem, I guess the next question is: can I use contours to create a DEM that will work better for this application?


texts:
All I'm trying to do is list out the feature classes in a geodatabase. As far as I understand if you use arcpy.ListFeatureClasses, it should return a list of the features classes in that geodatabase. But I get nothing. I don't get None. I don't get []. Simply nothing is returned and the module acts like it's ready for the next command. I've checked and there is a gdb under that name, in that place, that does have feature classes in it. What am I missing here?

import arcpy
arcpy.env.workspace = "L:\\gathr\\indonesia\\Sara\\VIIRS_Detections\\January2016.gdb"
arcpy.ListFeatureClasses()


texts: Is there a programmatic (or otherwise portable) way to convert E00 ArcInfo interchange files into ESRI format?


texts: why ee.Reducer.sum() and ee.Reducer.mean() get the same value when used in reduceResolution()

texts:
I'd like to reduce Landsat image resolution to MODIS. But when I use sum () and mean() reducer, the results are the same. Then I tried using Minmax, mode reducer, surprisingly, they all are the same. What's wrong? Here are my script and comments.
// those are from GEE
// Load a MODIS EVI image.
var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD13A1').first())
    .select('EVI');

// Display the EVI image near La Honda, California.
Map.setCenter(-122.3616, 37.5331, 12);
Map.addLayer(modis, {min: 2000, max: 5000}, 'MODIS EVI');

// Get information about the MODIS projection.
var modisProjection = modis.projection(); // what's the difference betwen .projection() and crs()
print('MODIS projection:', modisProjection);

// Load and display forest cover data at 30 meters resolution.
var forest = ee.Image('UMD/hansen/global_forest_change_2015')
    .select('treecover2000');

print('landsat projection',forest.projection());
Map.addLayer(forest, {max: 80}, 'forest cover 30 m');

// compare among different reducer///////////////////////
// Get the forest cover data at MODIS scale and projection.
var forestMean = forest
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(), // mean is 88.31
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });
// non of those works to compute class density
var forestSum = forest
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.sum(), // sum is 88
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });

Map.addLayer(forestSum, {}, 'forest cover sum at MODIS scale');

var forestMode = forest

    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mode(),
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });

Map.addLayer(forestMode, {}, 'forest cover mode at MODIS scale');

var forestMinMax = forest
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.minMax(),
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });

Map.addLayer(forestMode, {}, 'forest cover minmax at MODIS scale');


texts: How to convert a raster file (*.tif) into a composite geodataset by using Python in ArcGIS 10.2.1

texts:
I created geotiff rasters (approximately 6000 rasters) in .tif format by using Matlab Mapping Toolbox. When I examined the properties of the data, I realised that it is 64 bit double type. I need to create a batch process to extract values by using points and my code gave an error about the type of the data. (Error code 000860)

So it is obvious that I should use composite geodatasets in my process. Therefore I tried to convert them into composite geodatasets by using Composite Bands tool. When I used it for a sample data (which is opened in the table of contents) I was successful. But using this tool in a batch process gives an error which is:

Runtime error  Traceback (most recent call last):   File "", line 22, in    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 12855, in CompositeBands     raise e ExecuteError: ERROR 000271: Cannot open the input datasets Failed to execute (CompositeBands).

My code is given below:

import arcpy
import os
from arcpy import env
from arcpy.sa import *

mainfolder= "D:/csvfiles"
templateshp="D:/shp/station.shp"
mxd = arcpy.mapping.MapDocument("CURRENT")

for foldername1 in os.listdir(mainfolder):
    insidefolder=mainfolder+foldername1
    for foldername2 in os.listdir(insidefolder):
        insidefolder2=insidefolder+"/"+foldername2
        arcpy.env.workspace = insidefolder2
        for filename in os.listdir(insidefolder2):
            if filename.lower().endswith(".tif"):
                filename2=filename.replace(".csv.tif","")
                fullname=insidefolder2+"/"+filename
                fullout1=insidefolder2+"/"+filename2+".img"
                fullout2=insidefolder2+"/"+filename2+".shp"
                print "Raster Composite:" + fullout1
                arcpy.CompositeBands_management(fullname,fullout1+".img")
                print "Point conversion:" + fullout2
                ExtractValuesToPoints(templateshp,fullout1,fullout2,"NONE","VALUE_ONLY")
                print "OUTPUT: "+fullout
                for z in arcpy.mapping.ListLayers(mxd):
                    arcpy.mapping.RemoveLayer(df , z)
print "DONE!!!!"

What may be the problem is?

By the way, I am using ArcGIS 10.2.1.


texts: Help with Python script to query feature server

texts:
For the purpose of saving the data locally from a feature server online I have written code below:

import arcpy
fcout = r'C:\Forum\Hydrants\test01.shp'
url = 'http://gisweb.wsscwater.com/ArcGIS/rest/services/FireHydrant/FireBook/MapServer/0/query?text=&geometry=&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&objectIds=&where=1%3D1&time=&returnCountOnly=false&returnIdsOnly=false&returnGeometry=true&maxAllowableOffset=&outSR=&outFields=*&f=pjson'
fs = arcpy.FeatureSet()
fs.load(url)
arcpy.CopyFeatures_management(fs, fcout)

However only returns 1000 - there are 40000 data points, how would I create a shape file for the each of the other data points?


texts: Is the acronym GPS used by members of the public in countries like Russia or China, who have their own GNSS

texts: In the UK, and countries using the GNSS controlled by America, the acronym GPS (Global Positioning System)  is used to describe the data. I was wondering if members of the public in countries like Russia or China, who have their own GNSS, use a different term to describe GNSS data or whether, because of the proliferation of technologies using GPS, it is a term used across GNSS systems?


texts: How can I restrict the relation so the tables just join when is a relationship?

texts:
I made a relations (1:n) between 2 tables. Not every element from fathers layer have connection with elements from child layer. The connection is made between 'Tipo de Estrutura' Field and 'Est' Field like shown below.

The relation is made correctly. My problem is with the elements that have no relation (Example 'Rocha'). When I open the 'identify feature' of 'Rocha' Element is shown a relationship table with no elements. It is possible to only show relation when exist or hide this connection?


texts: Very slow connections to ArcSDE with over 20,000 datasets

texts:
I have a problem with ArcSDE that I have never come accross before and I was wondering if anyone could shed any light on this problem.

Problem ArcSDE 9.3.1:

  Connnecting to an ArcSDE instance (Oracle Spatial) with over 20,000
  datasets (raster) causes ArcGIS/ArcObjects code to connect very slowly - minutes in fact before a successful connection is made to the database.  However if a native SQL connection is made, the connection is almost instant.

Has anyone else had this problem and if so what can be done to remedy this problem.


texts:
I have a list of links and withing the links I have data attributes for lat and lon.

When I hover over the links I would like to show marker on the map.
let loc = document.querySelectorAll("a[data-lat]");

loc.forEach(node => {
    node.addEventListener("mouseover", locOver);
    node.addEventListener("mouseout", locOut);
})

function locOver(event) {
    var llt = event.target.dataset["lat"],
        lln = event.target.dataset["lon"];

    L.marker([llt, lln]).addTo(map);
}

function locOut(event) {
   //remove markers associated with the links
}

This seems to work partially but I need to figure out how to resolve a few issues:
Currently every time I hover over the same link it keeps adding more and more markers. I was thinking of removing the marker on mouseout event but later I found a function (setLatLng) that probably is more efficient and instead of adding and removing just uses the same marker to just change it's coordinates.
Also, when I am no longer hovering over links I'd like to remove the last link's marker. I do have a marker of my current position that I'd like to keep on the map.
I've tried it in many different ways but always get JS errors in console and all my solutions seem to be over-complicating the issue. I'm certain I'm pretty close, just stuck.


texts:
I am utilizing ArcGIS Pro and ArcPy to automate exporting of JPEGs from several counties in Texas. I am querying each county, and attempting to export each.
My script works fine when I run it in two separate part, but if I attempt to run them as one script only the basemap is exported, and not the polygons. See below.

import arcpy
import os

aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.listMaps()[0]
lyr = m.listLayers("myselection")[0]
lyt = aprx.listLayouts("Layout")[0]
mf = lyt.listElements("mapframe_element")[0]

county_names = []

with arcpy.da.SearchCursor('myselection','NAME') as cursor:
    for row in cursor:
        county_names.append(row[0])

    del cursor, row

for county in county_names:
    arcpy.SelectLayerByAttribute_management("myselection", "NEW_SELECTION", "NAME = " + "'" + county + "'")
    mf.zoomToAllLayers()
    mf.camera.scale *= 1.1
    arcpy.SelectLayerByAttribute_management("myselection", "CLEAR_SELECTION")
    lyr.definitionQuery = "NAME = " + "'" + county + "'"
    lyr.visible = True
    mf.exportToJPEG(path + county, resolution=300)
    lyr.definitionQuery = ""


texts:
I created a collection that contains a multipolygon document.

I want to write a query in spring data that tells me if a point (lon, lat) intersects the multipolygon.
The query in mongoDb is as follows:

db.my_multi_polygon.find({
    geometries : {
        $geoIntersects : {
            $geometry : {
                "type" : "Point",
                "coordinates" : [ 12.704157, 41.720075 ]
            }
        }
    }
});

but I would like to understand how to translate it into Spring Data

Is it possible to write a query that returns a boolean (true = intersects, false = not intersects)?


texts: gdalbuildvrt vs gdal_merge on overlapping images

texts:
I once have used gdalbuildvrt on georeferenced images. That was okay when images did not overlap, It was also time efficient compared to gdal_merge or gdal_warp. Now I have overlapping images and gdalbuildvrt is a little problematic to me. Let me explain:
For the testing purposes I have one image named 0.jpg(1600*900). I use the following command to georeference 0.jpg.

gdal_translate -a_srs EPSG:4326 -co interleave=pixel -a_ullr -180.0 90.0 180.0 -90.0 0.jpg ONE.tif

To see how it will merge overlapping areas I type the following:

gdal_translate -a_srs EPSG:4326 -co interleave=pixel -a_ullr -90.0 90.0 90.0 -90.0 0.jpg TWO.tif

So what I expect when I merge ONE.tif and TWO.tif is TWO.tif fills the center of the picture. Finally, I use gdalbuildvrt on these two:gdalbuildvrt -input_file_list LIST.txt ONETWO.vrt, and then gdal_translate ONETWO.vrt ONETWO.tif

My LIST.txt consists of two lines: first ONE.tif then TWO.tif.
The problem is here. 0.jpg was 1600*900 pixels. So I expect ONETWO.tif to be 1600*900 pixels-- ONE.tif covering from -180 to +180, and TWO.tif overwriting from -90 to +90. Although TWO.tif covers the center in the result, the final image is 2133*900 pixels. Why is this happening?  Of course,  I can use -outsize on gdal_translate to restore the original resolution but I am not sure if it is the right way to solve this.

On the other hand, gdal_merge is also an alternative, giving me what I want(not changing the original resolution) but it is likely that it will take longer compared to gdalbuildvrt.

I use GDAL 2.2.0, released on 2017/04/28 on OSGeo4W64.
Below is an example information stored in ONETWO.vrt. It is probably this area that needs to get modified.

<VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Red</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">ONE.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1600" RasterYSize="900" DataType="Byte" BlockXSize="1600" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="1600" ySize="900" />
      <DstRect xOff="0" yOff="0" xSize="2133.33333333333" ySize="900" />
    </SimpleSource>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">TWO.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1600" RasterYSize="900" DataType="Byte" BlockXSize="1600" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="1600" ySize="900" />
      <DstRect xOff="533.333333333333" yOff="0" xSize="1066.66666666667" ySize="900" />
    </SimpleSource>
  </VRTRasterBand>

Why are the destination coordinates messed up like 533.3333 or 2133.3333? And, is there a way to use gdalbuildvrt on overlapping images just like gdal_merge?


texts: Is it possible to create a GeoServer raster (GeoTIFF) data store where the URL to the actual file points to a web location instead of a folder in the GeoServer?


texts: ArcGIS Raster Mosaic Dataset and GDAL VRT, when to use?

texts:
ArcGIS Desktop has a the Raster Mosaic Dataset, QGIS and friends have Virtual Raster Tables.

A skim of the documentation reveals much similarity. They both create or enable an alternate representation of data that sits somewhere else without having to actually create a new image. Examples include creating a single seamless image from many tiles without actually creating that single image, or mix-and-matching channels or selectively removing value ranges from the result, and more.

How far does that similarity extend? Does an RMD enable things that VRT can't? What are they? What can one do with VRTs that you can't with RMDs?

In short, assuming one has equal access to and comfort using both tool chains, when should a person use Raster Mosaic Datasets and when are Virtual Raster Tables the order of the day?


texts: Spatial Join unique values for closest target feature

texts:
I have a polygon shapefile of land plots with their corresponding cadaster number and a point shapefile with a water meters for each lot.

What I'm trying to do is to join the cadaster number from the polygon shapefile to its corresponding water meter. For this I'm using a spatial join with the point file (water meters) as my target feature and the polygon file (land plots) as the join feature and CLOSEST as the match method.

The problem is that the water meters and the land plots do not allign correctly as shown in the figure below and it's adding the same cadaster number to two different water meters.

Is there a way I can tell the program to assign a cadaster number only once and if a second point has this same polygon as the closest one to move on to the second closest?

I have some experience with python so if there is a way to do  it this way it won't be a problem.

The drawing attached is just a simple example, what I actually have is this problem for a whole city so I can't just move the water meters manually.


texts: How to know the date when an image has been capturing on Google Hybrid QGIS?

texts: How can I know precisly the date of different Google Hybrid pictures used in my project in QGIS?


texts: ArcMap: error when selecting attribute records

texts:
At 10.8.1 I need to select all records from a shapefile string attribute field where the left-most character is blank (that is, a single whitespace).  For guidance, I rely on the onboard ESRI Help documentation by searching for "SQL reference for query expressions" (screenshot below). Following that Help example, I use the following expression:
LEFT (field_name, 1) = ' '

Unfortunately the selection fails with the error

An invalid SQL statement was used.

When I surround the field name with double-quotes (not two single quotes), as in:
LEFT ("field_name", 1) = ' '

The same error occurs.
This should be so simple!  Is the Esri documentation wrong? What am I missing?


texts: stack() in R does not include all .tiff files in my directory

texts:
I am trying to perform stack() in R using enhanced vegetation index (EVI) time-series maps but it doesn't include all the .tiff files in my directory. Here's the list of files when I dir(files.evi):

  [1] "MOD13Q1_EVI_2019_001.tif" "MOD13Q1_EVI_2019_017.tif"

  [3] "MOD13Q1_EVI_2019_033.tif" "MOD13Q1_EVI_2019_049.tif"

  [5] "MOD13Q1_EVI_2019_065.tif" "MOD13Q1_EVI_2019_081.tif"

  [7] "MOD13Q1_EVI_2019_097.tif" "MOD13Q1_EVI_2019_113.tif"

  [9] "MOD13Q1_EVI_2019_129.tif" "MOD13Q1_EVI_2019_145.tif"

  [11] "MOD13Q1_EVI_2019_161.tif" "MYD13Q1_EVI_2019_009.tif"

  [13] "MYD13Q1_EVI_2019_025.tif" "MYD13Q1_EVI_2019_041.tif"

  [15] "MYD13Q1_EVI_2019_057.tif" "MYD13Q1_EVI_2019_073.tif"

  [17] "MYD13Q1_EVI_2019_089.tif" "MYD13Q1_EVI_2019_105.tif"

  [19] "MYD13Q1_EVI_2019_121.tif" "MYD13Q1_EVI_2019_137.tif"

  [21] "MYD13Q1_EVI_2019_153.tif" "MYD13Q1_EVI_2019_169.tif"

However, when I do

evi <- list.files(files.evi, pattern = glob2rx("*MYD13Q1*.tif$"), full.names = TRUE)
all_evi_st <- stack(evi)

it only stack half of the files on my list:

  class      : RasterStack

  dimensions : 668, 1105, 738140, 11  (nrow, ncol, ncell, nlayers)

  resolution : 231.6564, 231.6564  (x, y)

  extent     : 11251086, 11507066, 1348703, 1503450  (xmin, xmax, ymin,
  ymax)

  crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
  +b=6371007.181 +units=m +no_defs

  names      : MYD13Q1_EVI_2019_009, MYD13Q1_EVI_2019_025,
  MYD13Q1_EVI_2019_041, MYD13Q1_EVI_2019_057, MYD13Q1_EVI_2019_073,
  MYD13Q1_EVI_2019_089, MYD13Q1_EVI_2019_105, MYD13Q1_EVI_2019_121,
  MYD13Q1_EVI_2019_137, MYD13Q1_EVI_2019_153, MYD13Q1_EVI_2019_169

Can anyone tell me what's happening and what's the workaround? I am new in remote sensing and still shabby with R.

Other info I can give about the time-series images is that they are from MODIS TERRA and AQUA time series derived from the 250 m 16-day composite vegetation index products (MOD13Q1/MYD13Q1).


texts: Can FGDC Status Checker be set to private?

texts:
I am currently using the FGDC Status checker to track the health and status of my ArcGIS Services.

https://statuschecker.fgdc.gov/

It works great, other than the fact that the test results and REST endpoint URLs are publicly available. Does anyone know if my account and test results can be set to private? Or of any other similar monitoring methods?


texts: Sorting of different size and color points (multiple attributes) from one shapefile in ArcGIS 10.3 possible?

texts:
I produced a map in ArcMap 10.3 with a point shapefile representing two types of values from its attribute table. One value is represented by symbol size, the other by color. For that I used the Symbology-->Multiple Attributes in the Layer properties dialog.
I can not find a way to sort the points, so that smaller points are in front and bigger in the background. Situation now is that some bigger points overlay smaller ones in areas with point clusters. The proposed methods I found online are only working for unique values, like sorting the data or using symbol levels.


texts: Filter Specific dates in Google Earth Engine?

texts:
I've called Landsat8 images from 2014 to 2019 using filter date function. But I don't know how can I call Landsat8 images only for June, July, and August from 2014 to 2019?

Map.centerObject(geometry);

var landsat = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterDate('2014-01-01','2019-01-01')
.filterBounds(geometry);

this is code's link:
https://code.earthengine.google.com/5c2959ec6adecfdefa817c37d74effe2


texts: QGIS 3 interprets polylines as multipolylines?

texts:
ogrinfo -al reports a file as Geometry: Line String and each feature as LINESTRING. But when loaded in QGIS 3, it is reported as Line (MultiLineString). As far as I remember, QGIS 2 imported it as LineString. Is this an intentional change or a bug?

(By the way, is the command "Multipart to Single Parts" in QGIS 3?)


texts: Identifying (and solving) polygon problem inside geometry in QGIS?

texts:
I've downloaded a countries shapefile from Natural Earth. In QGIS 2.14.11-Essen, in Vector|Geometry Tools|Check Geometry Validity... I got 2 errors (actually 1).

It says the error is in the "Feature 7". Opening the Attribute Table for this layer, I see numbers for each line, beginning from 0. The 7th line contains United Arab Emirates (UAE). So I suppose the problem lies in this country.

Giving a zoom to selection, I see this country in yellow, and notice it has eight small islands along its northern border. Also, it has a white patch on the extreme east of its territory (which belongs to Oman), and inside this patch, another patch, which again belongs to UAE (inside the red circle below).

First question: how do I know the number of the islands, to know which ones are the polygon 8 and polygon 9 given in the error message?

In this easy case, I just supposed the problem lied in that exclave, and after removing it (and the hole where it was contained), the problem was solved.

Second question: how do I solve the geometry problem, without removing that part of UAE's territory, and preferentially by hand, node by node (so I know exactly what's going on, what is not always the case when we use "automatic correction" tools)?


texts: No output in Python console when using PyQGIS

texts:
When I run my code in QGIS Python console I get output, but when I run my code from normal console I didn't get any output.
Why does it work only in QGIS Python console?
import sys, os

sys.path.append(r'C:/PROGRA~1/QGIS3~1.16/apps/qgis/python/plugins')
os.environ['PROJ_LIB'] = 'C:\\Program Files\\QGIS 3.16\\share\\proj'

from qgis.core import *
from qgis.utils import *

qgispath = r'C:/PROGRA~1/QGIS3~1.16/apps/qgis'

QgsApplication.setPrefixPath(qgispath, True)
qgs = QgsApplication([], False)
qgs.initQgis()

path = r'C:/Users/ja/Inne/Desktop/POLSKA/PRG_PunktyAdresowe_POLSKA.shp'
lyr = QgsVectorLayer (path, 'PRG_PunktyAdresowe_POLSKA', 'ogr')

for f in lyr.getFeatures():
    print(f.attributes())
    break


texts:
Continuously to:

SVG symbol colors differences in QGIS 2.12?

i tried to change the color of this pin mark

(i took the svg file from http://www.clker.com/cliparts/G/R/Y/M/S/g/map-pin-yellow.svg)

in point layer:

by opening the notepad of the svg file and adding those line into the notepad of the svg file

fill="param(fill) #FFF" stroke="param(outline) #000" stroke-
width="param(outline-width) 1"

but with no result.


texts: Solving coordinate differences between packages

texts:
I have two maps generated by two different methods. The first, I used maps package and the function map() and symbols(). On the other hand, I used the ggmap package with the ggmap() and geom_point() functions. This the structure of my data and my code:

# Capitalizing first letters
locations <- read.csv("ciudades.csv", header=TRUE, sep=";")
stri_trans_totitle(as.character(locations$ciudad))
write.csv(locations, file = "locations.csv")

# Obtaining latitude & longitude
lonlat <- geocode(as.character(locations$city))
ciudades_lonlat <- cbind(locations, lonlat)
write.csv(ciudadess_lonlat, file = "ciudades_lonlat.csv")

> str(city_frequency)
'data.frame':   144 obs. of  4 variables:
 $ city: Factor w/ 144 levels "A Coruna","Agen",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ freq: int  7 2 4 1 4 2 1 1 1 5 ...
 $ lon : num  -8.412 0.616 -8.446 -4.262 2.149 ...
 $ lat : num  43.4 44.2 40.6 42.8 43.9 ...

xlim <- c(-13.08, 8.68)
ylim <- c(34.87, 49.50)
map("world", lwd=0.05, xlim=xlim, ylim=ylim)
symbols(city_frequency$lon, city_frequency$lat, circles=sqrt(city_frequency$freq), bg="#93ceef", fg="#ffffff", lwd=0.5, inches=0.15, add=TRUE)

m <- get_map("France", zoom = 5, maptype = "toner")
ggmap(m) + geom_point(data = city_frequency, aes(lon = lon, lat = lat, size = sqrt(freq)), alpha = 0.7) + labs(x = "longitude", y = "latitude", title = "City Frequency")

The thing is that in the second map all the points are located a little bit North from the right place.

How can I solve this problem?

Note: in order to get the coordinates of my data frame I used the geocode() -which I have added at the first lines of code- in a early step.

---> I have created a public repository in github in order that everyone can access my data frames. Feel free to play with them.


texts: QGIS 2.8.6 raster calculator: cannot select cells which are below a value in a raster

texts:
I have a distance raster map in tif format and I am trying to generate a binary raster where cells below 200m are asigned value 1 and cells above that value are 0.
I've converted the distance raster to asc and checked that there are values under 200.
I've written the following expression in the raster calculator:
"disagua@1" <= 200
But the raster I get only displays 0s.
I've tryied to reclass the distance raster using Saga reclass tools but I just get a raster with 0 values.
However, If i try with values under 1000 I get a correct binary map.
Could somebody please help me?


texts: Setting bounds and making map bounce back if moved away

texts:
I want to make this map bounce to the center, and make the red line the limit that you can move the map.

Don't know why my fitbounds doesn't work.

Any help?

My code:

var mapBounds = L.latLngBounds([
                    [-8576 / 2, -8576 / 2],
                    [8576 / 2, 8576 / 2]
                ]);
                map.fitBounds(mapBounds);

http://plnkr.co/edit/XQxgpcYjiVaUhU4jCQV9?p=preview

I want something similar to this example:

https://www.mapbox.com/mapbox.js/example/v1.0.0/maxbounds/

You can't move away from that part.


texts: Building terrian data in ArcGIS Server for Cesium

texts:
I want to use ArcGIS terrain Data in Cesium. I don't know how to build terrain data in ArcGIS Server 10.3.1 according to specifications of Cesium.
Is there a tutorial or sample data which can be used in Cesium?
I know there is ArcGISTiledElevationTerrainProvider class in Cesium which can load ArcGIS terrain data but I don't know how we can build data for this provider in ArcGIS.


texts:
I use the gdal executable in QGIS file, with a DOS command.

I have a csv file like this (test.csv) :

  X       Y
-5.48   42.81
-4.78   42.52
-5.06   42.02

This coordinates are GCP. I would like to convert my csv file in a VRT. On the internet, I found example like :

<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>test.csv</SrcDataSource>
        <GeometryType>wkbPolygon</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <Field name="id" src="id" />
        <GeometryField encoding="WKT" field="geo" />
    </OGRVRTLayer>
</OGRVRTDataSource>

Simply, I don't know where to put this line. Should I put them in a DOS command in my QGIS file ?


texts:
I've been following this tutorial to server a very large raster tile set. The data set is many thousand 1km tiles covering the whole of London  (derived from LiDAR data) the data is stored on a PostgreSQL db the tiles are generated as the below command

raster2pgsql -s 27700 -t 100x100 -I -C -M {}/raster_*_3d_{}.asc -F solar_rasters.solarp_3d_{} | psql -d lsom

The error I receive when trying to add the store is

Could not list layers for this store, an error occurred retrieving them: Failed to create reader from file:data/mosaicpgraster.pgraster.xml and hints Hints: REPOSITORY = org.geoserver.catalog.CatalogRepository@40290072 EXECUTOR_SERVICE = java.util.concurrent.ThreadPoolExecutor@14db3d9f[Running, pool size = 0, active threads = 0, queued tasks = 0, completed tasks = 0] System defaults: FORCE_AXIS_ORDER_HONORING = http GRID_COVERAGE_FACTORY = GridCoverageFactory TILE_ENCODING = null COMPARISON_TOLERANCE = 1.0E-8 FORCE_LONGITUDE_FIRST_AXIS_ORDER = true FILTER_FACTORY = FilterFactoryImpl STYLE_FACTORY = StyleFactoryImpl FEATURE_FACTORY = org.geotools.feature.LenientFeatureFactoryImpl@45394b31 LENIENT_DATUM_SHIFT = true

In the log I can see;

20 Nov 15:50:28 ERROR [jdbc.custom] - No level available for oek
20 Nov 15:50:28 ERROR [imagemosaic.jdbc] - No level available for oek
20 Nov 15:50:28 WARN [imagemosaic.jdbc] - java.io.IOException: No level available for oek
java.io.IOException: java.io.IOException: No level available for oek

oek is present in the connection xml as shown in the tutorial, but I'm not clear what this value should be.


texts: Using FTP Connection with ArcGIS for Desktop?

texts:
I keep some data on an FTP server.

Can I conect with ArcMap 10.1 to use this data?


texts: No building and Point of interest available while setting up "Rendering Osm Data Ubuntu using MAPSERVER"

texts: I have followed RenderingOsmDataUbuntu 10.10 and set up Google styled web map using mapserver. The osm-google.map seems to be not using table osm_new_buildings. Also no POI (Point of interest) data is seems to be imported from the OSM file. How can I add both?