texts: Measure loading time of layers in Mapserver

texts:
I use mapserver to generate one single tile per request. The mapfile contains several layers and I would like to know which of them is the bottleneck that needs too much time to render.
WIth shp2img I can get a rough idea but as we also use mapscript it would be better to get these debug information directly from mapserver.
Is this possible?

shp2img -m alk.map -o ~/tmp/test.png -map_debug 3

[Tue Jun  4 11:13:34 2019].577359 msDrawRasterLayerLow(relief): entering.
[Tue Jun  4 11:13:34 2019].618554 msDrawMap(): Layer 0 (relief), 0.044s
[Tue Jun  4 11:13:34 2019].619952 msDrawRasterLayerLow(osm_ab_4000_sichtbar): entering.
[Tue Jun  4 11:13:34 2019].638350 msDrawMap(): Layer 2 (osm_ab_4000_sichtbar), 0.020s


texts: Python: ogr2osm create osm file from WFS

texts:
Does anyone know if ogr2osm can use WFS Layer as input source to create a osm file i could not find anything in the documentation expect that it should support all ogr file formats.

If so how can I do something like that (with ogr2osm)?


texts:
I'm trying to put an x/y offset to the point-geometry-layer to get map-tip-informations from all overlapping points like in the Point Displacement Renderer.

For getting map-tip-informations it seems to be important that the point is visible at the stored coordinate and is not modified by an offset in the renderer. So I have to recalculate the stored coordinates and modify it with a varying offset which depends for example on the actual scale. With this recalulated, displaced coordinates I'm trying to create a new layer. This displaced points should be recognized by map-tips.

I created a virtual fields (X_modified, y_modified) and virtual_geometry with the field name "geometry" (by "geom_from_wkt('Point(X_modified || Y_modified)') - but there is no geometry-field-type available.
The virtual layer didn't recognize the virtual geometry or other virtual fields...

-> Is it possible to create a virtual layer based upon virtual (geometry) fields of another layer ?

The calculation of the modified coordinates in the postgresql-data-provider seems not to be possible, because i do not know how to transfer the actual scale-variables (or other ones) from qgis to postgresql and its view-definitons. Or - do you know a way to do this ?
[QGIS 2.18]


texts: Spatial AND temporal join between two point data sets

texts:
I was hoping I could get some pointers in how to approach this issue, which is tricky version of a spatial join. I have been searching and I can't find a particular solution to what I need.

I have two data sets of points with related attributes that include date and coordinates. I want to join the points in data set A (records of events) with the nearest neighbor in the data set B (records of animal locations). The catch is that I want to find a neighbor in B that is not only closest to A in space but also in time. Perhaps explaining the purpose of this might help visualize what I want to do.

I want to join each event in A with the animal in B that was closest to the location of the event AND closest to the date in which the event happened. I have searched options with spatial join and "join by group" in ArcGIS. Yet I haven't been able to match according to attribute date given that the dates will not match exactly in every case. Thus, that is why I am looking for the record that is closest to that date instead of matching records that have equal dates and why I can't create fixed groups according to specific dates. Also I have a huge data set of about 100,000 records.

This is how my data bases look like. I have removed columns and rows to make it easier to visualize:

Database A, animal locations

indiv_code  sex month   year    Local_date  treatment   perturbati  xx  yy
A   M   3   2014    04-03-14    E   B   492533.82430000000  7935470.32100000000
A   M   2   2016    25-02-16    E   A   506713.28360000000  7938586.99800000000
A   M   2   2016    25-02-16    E   A   506732.55540000000  7938504.00800000000
B   M   2   2010    22-02-10    E   B   549766.00000000000  7909026.00000000000
B   M   2   2010    22-02-10    E   B   550720.00000000000  7908986.00000000000
B   M   2   2010    22-02-10    E   B   551746.00000000000  7906541.00000000000
B   M   1   2015    14-01-15    E   B   545292.81210000000  7891203.71800000000
B   M   1   2015    14-01-15    E   B   545292.82570000000  7891209.25100000000
B   M   1   2015    14-01-15    E   B   545262.86540000000  7891176.12900000000
C   M   4   2013    08-04-13    E   B   482280.00000000000  7941062.00000000000
C   M   4   2013    08-04-13    E   B   482266.00000000000  7941174.00000000000
C   M   4   2013    08-04-13    E   B   481805.00000000000  7942360.00000000000

Data base B, events

OID_    Dateincident    Year    Longitude   Latitude
1   10-01-08    2008    556517  7875861
2   15-01-08    2008    549887  7878550
3   15-01-08    2008    547800  7883200
4   17-01-08    2008    555500  7884500
5   30-01-08    2008    553530  7877180
6   31-01-08    2008    560328  7874984
7   31-01-08    2008    556517  7875861
8   01-02-08    2008    560400  7875900
9   01-02-08    2008    550400  7880500
10  01-02-08    2008    560500  7881500
11  01-02-08    2008    551800  7886800
13  14-02-08    2008    566300  7887700
14  20-02-08    2008    570326  7890088
15  27-03-08    2008    662271  7837444

I was thinking R could be the appropriate tool for this spatial and temporal join. Does anyone have any idea on how I could go about this? I have read so much information on spatial join, nearest neighbor and still can't figure this one out. If this is possible in ArcGIS, I have an advance licence for ArcGIS for Desktop.  If someone has an idea on how to do this in R and wants to help out, I can make a shorter sample of the data available to you.


texts: Is it possible to combine buffer zone analysis with time series analysis in QGIS?

texts: There is a multitemporal analyzer in QGIS where it is possible to do a buffer zone analysis combined with time series analysis, where there are two layers an input layer and a target layer and the effect of objects on the input layer can be analyzed the environment of the objects on the target layer specified distances? For example, we have a layer of a city road network with the data of car traffic in the roads, and there is a layer of roadside homes in the city with the data of cancer patients in homes, and can be analyzed that to air pollution on the road, on the input layer within what distances and to what extent it causes cancer on the target layer. How does the number of cancers change depending on the distance from the car's road? If we have carried out buffer zone analysis for several years, say in 2002, 2003, 2004 and 2005 then in the multitemporal analyzer these layers and the analyzed data of the years could be grouped, and time series analysis could be done about it that how the impact of road-side air pollution on cancers has changed at specified distances from year to year. Data should not only be analyzed across the whole layers but also within an Optionally selected part of the layers.


texts:
I manage a fairly large estate of QGIS installs in the UK and am looking to set the default date formats to the uk dd/mm/yyyy format from our MSSQL data warehouse.  Its set and works correctly in the database but QGIS defaults to YYYY-MM-DDT00:00:00.  Is there an option anywhere for me to change the format of this display date globally?

I am aware I can do it via widgets however we have over 1000 datasets in the warehouse and any local datasets the users have would still be affected by the default date formatting.


texts:
I have two layers, one is a point layer of schools and another is a polyline layer of streets. The end goal is to create a buffer along the lines that is 20 feet from the centerline that goes 1 mile in distance from the school/point layer.

Perhaps I am not phrasing this question correctly so if it is already answered somewhere else, I would gladly accept that.

I am using QGIS, which I am brand new to using. I have several years experience with ArcGIS (student license expired since graduation).

Would the plugin for QChainage accomplish this? Has anyone used this plugin?

If I had access to ArcGIS, it looks like this would be a Linear Referencing question. Would Locate Features Along Route be able to solve this problem?

Edit: Measuring the distance between lines and points in QGIS looks like it's finding the distance from a point to a line which is not what I am asking. I am asking to find the distance along a line from a point.


texts: Calling the correct element of an output via PyQGIS Script

texts:
I am trying to apply a style to the output file of a script that I´m running from the toolbox in QGIS 3.4.7. The last step is as described below.

        def processAlgorithm(self, parameters, context, feedback):

        if parameters['Input1'] is not None:
            source = self.parameterAsSource(parameters,'Input1',context)
            (sink, dest_id3) = self.parameterAsSink(parameters,'Output1',context,source.fields(),source.wkbType(),source.sourceCrs())
            total = 100.0 / source.featureCount() if source.featureCount() else 0
            features = source.getFeatures()
            for current, feature in enumerate(features):
                sink.addFeature(feature, QgsFeatureSink.FastInsert)
                feedback.setProgress(int(current * total))

After that, I would like to apply a fixed style to the output file so that it will be displayed on Layers with the applied style. I reached the lines below:

                self.source.loadNamedStyle(u"C:\\Users\\ct279359\\Desktop\\Styles\\Roads.qml")
                self.source.triggerRepaint()

I got this error:

  AttributeError: 'InOut' object has no attribute 'source'

I would like to know the correct way to call the output file created in the first part.

I did some research on the subject Here, Here, Here, Here, Here and Here.

All explain ways to invoke the QML file or how to change the style of a file that is already in the Layers section. I need next to this, to invoke the correct output file.

-------------------------------------------------------------
Full code here.

from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsFeatureSink,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterVectorLayer,
                       QgsProcessingParameterFeatureSink,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterRasterDestination,
                       QgsVectorLayer,
                       QgsProject)
import processing

class InOut(QgsProcessingAlgorithm):

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return InOut()

    def name(self):
        return 'InOut2'

    def displayName(self):
        return self.tr('InOut2')

    def group(self):
        return self.tr('Auxiliares')

    def groupId(self):
        return 'auxiliares1'

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('Input1','Input1',optional=True, types=[QgsProcessing.TypeVector], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Output1','Output1', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, feedback):

        if parameters['Input1'] is not None:
            source = self.parameterAsSource(parameters,'Input1',context)
            (sink, dest_id3) = self.parameterAsSink(parameters,'Output1',context,source.fields(),source.wkbType(),source.sourceCrs())
            total = 100.0 / source.featureCount() if source.featureCount() else 0
            features = source.getFeatures()
            for current, feature in enumerate(features):
                sink.addFeature(feature, QgsFeatureSink.FastInsert)
                feedback.setProgress(int(current * total))

                self.source.loadNamedStyle(u"C:\\Users\\ct279359\\Desktop\\Styles\\Roads.qml")
                self.source.triggerRepaint()

        return results


texts: Using GDAL Virtual I/O (/vsicurl) when Adding a Geoserver Store

texts:
I was under the impression that configuring Geoserver with the GDAL plugin would allow me to utilize GDAL's virtual I/O prefix for accessing data stores in a zip file, s3, etc. using /vsicurl/ or /vsis3/.

Having gone through the installation of GDAL and geoserver, I am able to use gdalinfo with the vsi prefix to successfully retrieve information about a file hosted in cloud, but when attempting to do this by adding a data store, I'm pretty sure Geoserver thinks this is just a file.

Any tips?


texts: How to work with 20 m spatial resolution bands in Sentinel-2 Imagery Level2A GEE

texts:
I am using the following script to make a query of Sentinel-2 Imagery Level-2 with GEE code editor.

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Finding images in the Sentinel2 L2A Collection
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
var s2_2018_03_01TO2018_05_31 = ee.ImageCollection("COPERNICUS/S2_SR")
  .filterMetadata('MGRS_TILE', 'equals', '32UQC')
  .filterDate('2018-03-01', '2018-05-31')
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 30));

print(s2_2018_03_01TO2018_05_31)

I attached an image showing the bands found in every image element of my image collection. I used SNAP in the pass to process the Level-1C Sentinel-2 Imagery to Level-2A. I can select which resolution I would like to work with. In my project I would like to use 20m spatial resolution imagery. Question, is there a way to have all bands in this resolution, instead that having elements with different spatial resolution in my image collection? Is there a method to scale down the bands B02, B03, B04 and B08 to 20m spatial resolution?


texts: Extracting number and attributes of points in a GRID in ArcGIS

texts: I will be thankful if someone can help in extracting the number of points with attributes using the grids. For example I want to know number of species points in a grid with attributes. I have ArcGIS and QGIS software


texts: Identify and remove isolated fishnet polygons

texts:
What process can I use to remove isolated cells in my layer; i.e. cells that don't
have an immediate neighbour are to be removed - I've circled a couple in the image below. These are based on an index within a table built with the fishnet tool. Many of the cells for the index column in the attribute table are blank.


texts: GeoServer - How to create source image from directory (JP2) and not only file by file

texts:
I would like to import all JP2 file from a directory with GeoServer.
I don't know how to do that. I can create image source with only one file JP2 but not with all file.


texts: How to save imagery from online raster service for offline use?

texts:
QGIS allows me to connect to free map services like ESRI's 'world imagery'.

I would like to save a certain part of it (decided by map extent) in the highest zoom level available in the service for later offline use.

this is for a purely non-commercial personal project. I intend to create a raster service for a few friends and myself in a connected environment. on this we want to overlay data related to some public service related work we are doing.

could someone tell me
a) if I can save the relevant raster images using QGIS
b) if yes, how to do it and in what format should I save those so that I can create a raster service from it


texts:
Does anybody have any experience with using Mapmanager to manage their MapServer mapfiles? My question is specifically about using INCLUDES in mapfiles. Since having an INCLUDE in a mapfile results in error message saying "Unable to access file. Error opening included file". Which leads me to believe it is supported and I am probably just doing something wrong.
Edit:
Tried using both relative and absolute path. The way I have my mapfiles is that I have them split multiple separate files MAP - contains the main map tags, LAYER - contains all the data, styles etc., CONFIG - contains any other potential information like db connection. MAP reads LAYER & CONFIG files using INCLUDE.  By themselves it works. But when I import the MAP file into mapmanager I get the error message. Any advice or tip on how to get around this would be greatly appreciated.


texts: ArcGIS Data driven pages Two MXD in one pdf

texts:
I have two mxd files created by data driven pages.

How is it possible to export them to one pdf file?

First a page from the first mxd file, then the first page of the second mxd file.
Next the second page from the first mxd file, then the second page of mxd file.


texts: How do I create a quick test/demo app with ArcGIS Mobile

texts: We are trying to do a quick evaluation of ArcGIS Mobile 10 on Windows.  Part of that is determining whether we can load online maps (like arcgisserveronline for example) and what the performance/rendering quality is like.  We tried using the 9.3 TileServices example but it does not appear to work with 10.0.  Anyone else tried anything like this?  Anyone know any good links on how to quickly get started on displaying maps on mobile (without our own map server)?


texts: Dealing with no-data values in raster calculator of QGIS

texts:
I am trying to patch two rasters with raster calculator but I am not finding the way to work with no-data values.

Lager@1 has some information I need as first option, this has as no-data -9999 (reading in the properties)
Lager@2 has the information I want to complete what is nodata in Lager@1

This command does not produce values where Lager@1 is no-data

("Large@1" != -9999) * "Large@1" + ("Large@1" = -9999) * "Lager@2"

I would like to avoid or skip the reclasification becase I might have many zeros and I need to perform this with diferrent rasters that might have different no-data formats, for instance Lager@2 has -99999 instead of -9999.

In GRASS this would be an r.patch task, but I cannot (or I do not know how to) control the order in the r.patch available in the Processing Toolbox.

r.patch input= Lager@1, Lager@2 output=Lager@3


texts: Remove a layer from TOC in my printed PDF?

texts:
I am new at writing Python code so I'm just grasping the basics.
I am trying to clip two shapefiles so that they are within the map boundaries. After that, I want to remove the unclipped layers and add the clipped layers to the TOC. Finally, print the map to a PDF.
As of now, my clipping seems to work, but when the PDF is printed, the layers I do not want on the map are still visible.

import arcpy
from arcpy import env

#Local Variables

RI_Schools = "E:\\Python Lab 2\\PythonInArcGIS\\RhodeIsland.gdb\\RI_Schools"
RI_Sewers = "E:\\Python Lab 2\\PythonInArcGIS\\RhodeIsland.gdb\\RI_Sewers"
Providence = "E:\\Python Lab 2\\PythonInArcGIS\\RhodeIsland.gdb\\Providence"
Schools_shp = "E:\\Python Lab 2\\PythonInArcGIS\\Schools.shp"
Sewers_shp = "E:\\Python Lab 2\\PythonInArcGIS\\Sewers.shp"

#Process: Clip

arcpy.Clip_analysis(RI_Schools, Providence, Schools_shp, "")
arcpy.Clip_analysis(RI_Sewers, Providence, Sewers_shp, "")

mxd = arcpy.mapping.MapDocument("E:\\Python Lab 2\\PythonInArcGIS\\Providence.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]

newLayer1 = arcpy.mapping.Layer("E:\\Python Lab 2\\PythonInArcGIS\\Schools.shp")
print newLayer1
arcpy.mapping.AddLayer(df, newLayer1)

newLayer2 = arcpy.mapping.Layer("E:\\Python Lab 2\\PythonInArcGIS\\Sewers.shp")
print newLayer2
arcpy.mapping.AddLayer(df, newLayer2)

arcpy.RefreshActiveView()
arcpy.RefreshTOC()

#Remove RI_Schools and RI_Sewers from TOC

names = [RI_Schools, RI_Sewers]

df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
layers = arcpy.mapping.ListLayers(mxd, "*", df)

for layer in layers:
  if layer.name in names:
    layer.visible = False

arcpy.RefreshTOC()
arcpy.RefreshActiveView()

#Export to PDF

arcpy.mapping.ExportToPDF(mxd, ""E:\\Python Lab 2\\PythonInArcGIS\\Providence1.pdf")
del mxd


texts: Watermarking photos using geotag and creating shapefile from coordinates?

texts:
I have geotagged photographs and I would like to place a watermark of the photos location on the front for easy reference as well as creating a shapefile from the coordinates.

Does anyone know of any open source software which could do this?

This link is an example of what Iam trying to acheive


texts: ArcMap ModelBuilder has tools blank on editing model?

texts:
When opening previously made models the tools appear blank. I cannot enter/edit the tools.

I run on ArcGIS 10.4 with a basic license.

The installation is new and the problem is persistent with all models.


texts: Calculating distance and/or time taken to travel from population centre to nearest hospital in QGIS

texts:
I'm trying to find out how long it would take to travel by car from each of the Lower Layer Super Output Areas (LSOA) in Norfolk, England to the nearest Hospital (or just what the distance is and I can estimate drive time, if the former is too difficult/impossible!).

Any suggestions how I can do that using QGIS, OSM road data and/or Ordnance Survey Meridian 2 road data, population weighted centroids for each of the 2011 Census LSOA areas and hospital locations?

I have all of the above in QGIS (see image below).

There are 537 LSOA and four hospitals in total. Is there a way to work out the time taken and/or distance to travel without having to do it one by one?


texts: Clipping feature layer to polygon to find total population that resides within polygon area using ArcGIS Desktop?

texts:
I am using ArcGIS-10.3.1.  I have a feature layer that includes total population by census tract (and includes lots of other info like shape length, shape area, population density, etc.).  I have an 11-km square polygon drawn on the map (just using the ArcMap drawing tool).  I want to know the total population that resides within the 11-km square.  The square includes some whole census tracts within it, but cuts across some others.

I assumed I could use the Clip tool somehow, but I'm unsure if it's even possible to clip a feature layer based on a polygon.  And, I'm unsure how to turn the polygon into a feature layer to then use it with the Clip tool.

Any ideas on how I could do this?  I'm a complete beginner - I can make maps, but don't know how to use Tools well.


texts: Hiding layers in mapfile (Mapserver) from WMS requests

texts:
We have a Mapserver mapfile used in a web app with a lot of layers. We want to publish only some of  them with WMS, keeping the others hidden to any WMS request.

It has to be the same mapfile and we need each layer has its NAME attribute specified

Is there any way to do this?


texts:
I have a polyline layer of roads and a raster of land ownership. For each raster cell I would like to calculate the distance to the nearest road.

I am on ArcMap 10.2 with Spatial Analyst and 3D analyst.


texts:
I am using QGIS 2.14.0 and I need to fix some topology errors. Problem is that when I move polygon node to the another node they seem to snap correctly. However, when I move that node again, in too many cases snapped nodes don't move together, and of course topology testing will fail. Is there any way around this problem?

I attached few images below. I've tried to turn off and on "Enable snapping on intersection" and "Avoid intersections". I've also experimented with different snapping tolerances with no success.


texts:
I'm trying to embed a PDF document into a Leaflet map popup, so that the PDF displays within the popup balloon when the user clicks the map point.

I've attempted to place <embed src="mypdf.pdf" width="500px" height="500px" /> into the leaflet popup wrapper elements, but the <embed> element does not show the content of the PDF, it is just blank.

However, I can take this same <embed> element and place it anywhere else on the web page, outside of the Leaflet popup element, and it shows the PDF file as expected - pointing to a problem specific to the popup.

Is there a reason why the <embed> tag is not working inside of a popup?


texts:
I'm trying to do Dark Object Subtraction (DOS) for a Spot2 and 4. I have calculated the TOA reflectance  and now want to get the at surface reflectance. Using QGIS

I would like to know how I can do DOS in QGIS, can I build a model or use an equation to do this. I have tried doing 6S however SPOT2 and SPOT4 HRV1/2 is not listed in the Pre-defined satellite bands.


texts: I am quite new to GIS software. I use QGIS. I have a question that might be trivial. I would like to derive the coordinates (Y and X) for an area. Just to make an example: I have two regions A and B. From the shapefile I know which are their coordinates. However, I would need the coordinates for the unique region A+B. Could you please give me some hints on how to do that?


texts: How to project TIGER & WGS84 to Spherical Mercator?

texts:
So I've been banging my head on this for about a week now. I'm using QGIS on a Mac. I have various shapefiles in WGS84 and Census TIGER NAD83 (EPSG:4269) and I would like to give them the same projection as Google, Bing and Openstreetmaps.

My map of the world, the one in WGS84, is always slightly off. Here is a screenshot of Florida compared to a Cloudmade tileset:

The map on the bottom is Cloudmade and the grey one on top is mine.

I'm starting to think that I could be doing this wrong in QGIS. Here are my steps:

Import shapefile.
Settings -> Custom CRS -> Create & save a custom CRS called Spherical Mercator and using parameters: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
Set the layer and project to be using that CRS.
Right clicking on layer and save as shapefile.

I'm using Tilestache to display my tiles. They are getting rendered through Mapnik. The Map XML property "srs" is also set to: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

Any insight on what my issue could be would be greatly appreciated.

Thanks.


texts:
I am trying to add a converted shapefile (now a geojson) to a Leaflet basemap, using EPSG 26986, so I tried using proj4leaflet to reproject it. However, I keep getting an error message "phi2z has NoConvergence" - in fact, I seem to be getting one message for each point feature in the shapefile (there are over 800). I'm not sure why I'm getting this message - anyone have any thoughts?

Proj4js.defs["EPSG:26986"] = "+proj=lcc +lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs";
var map = L.map('map').setView([42.317007, -71.088560], 12);
L.tileLayer('http://{s}.tile.cloudmade.com/655da89c06dd43a7aee306e594341325/997/256/{z}/{x}/{y}.png', {
    attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © CloudMade',
    maxZoom: 16
}).addTo(map);

$.getJSON("image_locations.geojson", function(data) {
    data.crs={
        "type": "name",
        "properties": {"name": "urn:ogc:def:crs:EPSG::26986"}
        };
    L.Proj.geoJson(data).addTo(map);
});

live version here: http://altani.continuation.org/trilobites/leaflet.html


texts: I get the above error when trying to add a new raster or vector layer. Does anyone know how to correct it as the resulting projection is nowhere near where I need it.


texts: Turn all fields on in layer properties using ArcPy?

texts:
I've been trying to turn all fields of a layer through properties fields

I just need to turn them all on

 import arcpy
 from arcpy import env
 import sys
 import os
 mxd = arcpy.mapping.MapDocument("CURRENT")
 df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
 for lyr in df:
     if lyr.name == "roads":
         lyr.visible = True
         arcpy.RefreshTOC()
         arcpy.RefreshActiveView()
##Then here I want to turn all fields on of the layer name roads


texts: Resources for Calculating the Geometry In Oblique Photographs

texts:
Currently I am using supervisory classification algorithms to calculate the landuse within a certain image, the issue being the images are sourced from high oblique photography. I know the focal length and dimensions of the camera used, and the GPS position of where the photographs where taken.

I am looking for resources that will enable me to calculate the areas of the landuse from the information I already have. The idea being that if knowing the geometry of the photograph, I will be able to work out the relative distance between each pixel. Thoughts?


texts: Insert an input variable in a query in PostgreSQL

texts:
I'd like to insert an input variable (integer) in a query SQL (@variable) using the algorithm postgresql execute sql in the graphical modeler of QGIS.
For example, in this kind of query
INSERT INTO table2 SELECT * FROM table1 WHERE ID = @variable;

where I should insert the variable in the input panel of the modeler and click play to do the query.
I used pre-calculated value but with no success. Any suggestion?


texts: Calculating population within a polygon from a raster layer

texts: I am trying to calculate the population within a set of polygons (representing districts). I have a raster file of population data. Can I assign the population data to the polygons and then sum the population?


texts: Integrating User Defined Coordinate System into PyQGIS?

texts:
Here is my current code:

#Loading of post-wrf(raster) files starts here
crs = 'USER:100000'
QSettings().setValue('/Projections/defaultBehaviour', 'useGlobal')
QSettings().setValue('/Projections/layerDefaultCrs', crs)
rLayer = QgsRasterLayer("/home/francis/QGIS/tiffout/test.tif", "")
if rLayer.isValid() is True:
    print("Layer {} loaded successfully".format(rLayer.name()))
else:
    print("Layer {} did not load successfully".format(rLayer.name()))
QgsMapLayerRegistry.instance().addMapLayers([rLayer])

#Loading of Shp file (vector) starts here
crs = 'EPSG:4326'
QSettings().setValue('/Projections/defaultBehaviour', 'useGlobal')
QSettings().setValue('/Projections/layerDefaultCrs', crs)
vLayer = QgsVectorLayer("/home/francis/QGIS/shpfile/PHL_adm1.shp", "Map", "ogr")
if vLayer.isValid() is True:
    print("Layer {} loaded successfully".format(vLayer.name()))
else:
    print("Layer {} did not load".format(vLayer.name()))
QgsMapLayerRegistry.instance().addMapLayers([vLayer])

my .shp file has the default coordinates of EPSG:4326
my raster file has a user defined coordinate and I can't save these two layers as one image.

I am new to using QGIS and Python.

The first screenshot is the sample output that I want to achieve with the correct alignment of the coordinates. (The layers were manually loaded in the QGIS UI)

The second screenshot of the User Defined Coordinate that QGIS provided. and my .gr2 data is an output of WRF


texts: Complex 3D/elevation profiles in ArcMap 10.6

texts:
I'm looking to create a more informative 3D Profile than appears possible with vanilla Arc tools. Honestly I'm kinda surprised that this doesn't exist, maybe I'm searching the wrong terms.

Situation:
I have a Pipeline made up of sections, there are many sections, topology is touching, but it's not a network. At intervals along the pipeline there are fittings (stored as 3D points, snapped to the polylines).

Outcome:
I'd like a image showing the elevation profile of the pipeline selected, the fittings along the pipeline, any intersecting polylines and various attributes (ID, elevation, chainage).

Something like this:

Currently I'm doing this manually using GlobalMapper (which displays lines crossing the input on it's 3D profile) and PowerPoint (labels).

I'd like a way to do this within Arc, or with python, in such a way it can be incorporated into data driven pages, or at least automated to the point I can give a script a whole network and have it dump out sanely-named jpegs.

Environment:

ArcMap 10.6, basic and standard licences, 3D Analyst, Spatial Analyst
1m DEM
3D polylines and points in a FGDB
Reasonably strict IT system - additional python libraries are okay, a full software package is okay, adding VBA to Arc is a PITA

Things I've tried:

Asking ESRI
Searching here and on Google
Searching ESRI's community pages
QGIS plug-ins

I found:

Xacto Cross Section Tool (a possibility, but requires VB libraries)
https://crossviewgis.com/ (a possibility, but I have to exhaust other avenues first)
https://geogeek.xyz/how-to-make-a-fancy-elevation-profile-including-xy-points.html - very, very close to what I want, but as manual as my current workflow

I'm not scared of python, but I'll need a little more than the name of a library and some pseudo code to get a result.

We have QGIS, again I'd need a reasonably detailed run down of ant proposed process.

I wondered if there was a way to fake things by adding an XY plot to the profile?

Anything considered, even "that's too hard, don't bother"


texts: Tiled Web Map; what tile provider offer antarctic map tiles?

texts:
General Idea: I plot Sentinel-2 and Sentinel-3 images on a web map in a Jupyter Notebook using ipyleaflet. This works perfectly fine with the default Mercator projection. For images at high latitudes I use the Arctic Web Map as follows:

tls = ipyleaflet.TileLayer(opacity=1.0,
                           url='https://{s}.tiles.arcticconnect.org/osm_3575/{z}/{x}/{y}.png',
                           zoom=0, max_zoom=10)

m = ipyleaflet.Map(default_tiles = tls, zoom=6, center=[0,0])

This lead to a good result. Now I would like to do the same for the South Pole. Does anybody know a tile provider for projections like EPSG:3031?


texts: Getting started with MapServer 6 and Google Maps

texts:
I'm relatively new to GIS, but what I'd like to do is display map layers on top of Google Maps.

From what I understand, MapServer provides the tiles that can be overlaid onto a Google Maps page.

I'm proficient with HTML, JavaScript, Unix, and PHP so technically this shouldn't be a problem, but I'm finding a severe lack of resources when I am trying to find out how one gets started. The MapServer documentation seems to be a little over the place, and Google results aren't helping.

Where does one get started with GIS and mapserver with a view to integrate on top of Google Maps?


texts: How to assign the same point to duplicate records in ArcMap geocoder?

texts:
What I am trying to accomplish is to geocode a long list of addresses by assigning the point manually. The problem is that there are many duplicates which can not be removed. Therefore, I need to assign the same point to several duplicates causing the overlapping problem!

A quick discussion with my colleagues brought to my attention that this problem is an old problem non solved by any updates from ESRI. I am using ArcMap 10.2 and I was wondering if there is any way to get a shortcut on this problem even with a python script or else.


texts:
When I run the maximum likelihood classification it shows that the task has been performed but there is no output. However, 3 files can be seen in recycle bin namely input, output and error every time I run the tool.

There is no error message seen.

I'm not saving the output starting with a number as name. The name that I'm using has less than 8 letters

Software: ArcMap; Version: 10.6


texts: Area of Intersection between two polygon files in QGIS

texts:
I'm new to mapping concepts and software.

I need to report on the intersecting area of postcode polygons and NBN Rollout areas. So postcode 3000 has 3 units of area in rollout region 1 and 6 units of area in rollout region 2.

Postcode polygons come in an ESRI Shapefile.
NBN Rollout regions come in a KML file.

I've loaded both into QGIS, saved them into a new file with CRS = EPSG 4326 then loaded these new files. Now the layers have the same CRS.

I then went to Vector > Data Management Tools > Join attributes by location.

Here are the settings I've selected:

That runs for a few minutes before throwing this error:

The log messages panel gives the following log entry:

  Uncaught error while executing algorithm          Traceback (most recent
  call last):
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py",
  line 203, in execute
                self.processAlgorithm(progress)
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\algs\qgis\SpatialJoin.py",
  line 158, in processAlgorithm outFeat.setGeometry(f.geometry())
  TypeError: QgsFeature.setGeometry(QgsGeometry): argument 1 has
  unexpected type 'NoneType'

The final output I want to see is something like this:

postcode    rollout_region  area
3000        1               27
3000        2               44
3141        4               94
3141        6               46
3006        8               70
3006        9               41
3006        11              38

if the problem can be solved more easily in some other open source software I'm happy to try that.


texts: Stitching lines / endpoints with postgres within a certain radius

texts:
I have lines from an OpenStreetMap database, which I selected like this:

SELECT sq2.osm_id, sq2.type, sq2.simple_geom, sq2.z_order FROM (
    SELECT sq1.osm_id, sq1.type, sq1.z_order, ST_UnaryUnion(sq1.simple_geom) as simple_geom FROM (
        SELECT osm_id, type, z_order, ST_SimplifyPreserveTopology(geometry, 50.0) as simple_geom FROM osm_roads_gen0 WHERE "type" = 'primary') as sq1
    GROUP BY osm_id, type, simple_geom, z_order) as sq2
    WHERE ST_Length2D(sq2.simple_geom) > 150.0
GROUP BY osm_id, type, simple_geom, z_order;

This gives me the result below - lines where I have a bunch of lines, but the highway departures are discarded (lines below 150m, I don't need them). Now what I wanted to do is to "stitch" these lines (by their ref). Meaning I want to either modify their endpoints or insert a line between and use ST_LineMerge.

How can I do this using PostGIS? Something like ST_Merge_Endpoints(ST_DWithin(geometry[i], geometry[i + 1], 150.0))

The second problem would be to move the line ending points (after the merge) to the nearest line (to "reconstruct" the simplified highway departures). I can get the ST_Endpoint, but then what? How can I move the endpoint to the nearest line?


texts: Difficulty in getting QGIS to connect with in-built GPS (COM3 NMEA port)

texts:
I have a Motion J3500 rugged tablet which was upgraded last year to run Windows 10.

The tablet has a built-in GPS / broadband module (Gobi Sierra 2000).
I use QGIS (v3.2) on the tablet, with the GPS plugin installed.

I want to use the GPS to add data to my QGIS by using live tracking while out on site.

The GPS itself functions fine (see photo). When I open the GPS pane in QGIS and click connect (using 'find settings automatically'), it appears to have connected, but the GPS pane within QGIS is not showing any satellite and positioning data - just zeros.
The GPS's own window shows that the GPS is working (see photo), but this is not reflected in the corresponding pane within QGIS.

One thing I've noticed is that the GPS is outputting to the COM3 NMEA port, rather than COM5 or COM6 - and I seem unable to change this.

It could be a settings or drivers issue somewhere, does anybody have experience of this and how to fix / work around?

PS, I have used the GPS to GIS link successfully on this machine in the past, but that was before the Windows 10 upgrade from W7, and was with ArcGIS, not QGIS


texts: ArcPy save edits without stopping an edit session

texts:
I am editing a versioned SDE Feature class with arcpy (ArcGIS 10.2.1).

I know that to edit versioned data you should use the arcpy.da.Editor, start editing, start operation, do your edits, stop operation, and stop editing (use stopEditing(True) to stop editing and save edits).

My question, is there a way with arcpy to save your edits without stopping your edit session like you can in ArcMap?

If you right click on "Editor" then "Save Edits" on the Editor toolbar in ArcMap it will save your edits and not stop your edit session.

Is this possible using arcpy?


texts:
I have downloaded the following precipitation data in .nc4 format(attached). I have opened the 'PRECIPITATION' variable in the file in QGIS and I have made a note of its statistics. As working with nc files consumes more data and time, I have written the data of precipitation variable in an ascii file using MATLAB. I have opened the ascii file in QGIS and I have found that the statistics i get from the properties->metadata is different from that of nc file and if i run the zonal statistics of the ascii file for the entire polygon boundary, the statistics are matching with the ones of nc file. Can anybody please help me with the reason behind this.

ASCII File
https://drive.google.com/file/d/0B8OJc_3nbFeLcTlWbTE1YzA4TUU/view?usp=sharing

nc file https://drive.google.com/file/d/0B8OJc_3nbFeLTDM3Z0lXQktUVGFxaTl6U1NZc3EzcFhkUzRR/view?usp=sharing


texts: How to overlay subset of one data set onto another

texts:
I have two data sets: one for Thurston County, Washington and one for all of Washington state. I created a map of hydrology features in Thurston County, and now I want to take the data from the Washington state data set and overlay it on top of the Thurston County data set, so that all of the WA state data points inside the county boundary are added, but nothing else is. How would I do this in ArcGIS for Desktop?

I'd be interested in instructions on how to do this, but just knowing if there is a name for this type of operation so I can search for more references would be helpful.


texts: SRTM HGT to GeoTIFF conversion using GDAL

texts:
I have (for example) SRTM-1 file n48e083.hgt with resolution 3601x3601.
What is right way to convert it to GeoTIFF?
gdal_translate -srcwin 0 0 3600 3600 -a_ullr 83 49 84 48 n48e083.hgt n48e083.tif

or
gdalwarp -te 83 48 84 49 -ts 3600 3600 n48e083.hgt n48e083.tif

I suppose the first is right and faster but I have some doubts.
And second question.
I wonder why gdal_translate does not copy "Unit Type" metadata value from hgt to tiff but gdal_merge does in the examples above?
UPDATE 1: gdalinfo n48e083.hgt
Driver: SRTMHGT/SRTMHGT File Format
Files: n48e083.hgt
Size is 3601, 3601
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (82.999861111111116,49.000138888888891)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (  82.9998611,  49.0001389) ( 82d59'59.50"E, 49d 0' 0.50"N)
Lower Left  (  82.9998611,  47.9998611) ( 82d59'59.50"E, 47d59'59.50"N)
Upper Right (  84.0001389,  49.0001389) ( 84d 0' 0.50"E, 49d 0' 0.50"N)
Lower Right (  84.0001389,  47.9998611) ( 84d 0' 0.50"E, 47d59'59.50"N)
Center      (  83.5000000,  48.5000000) ( 83d30' 0.00"E, 48d30' 0.00"N)
Band 1 Block=3601x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768
  Unit Type: m

UPDATE 2:
gdalwarp with "-r cubic" option:
gdalwarp -te 83 48 84 49 -ts 3600 3600 -r cubic n48e083.hgt n48e083.tif

makes tiff file which appears to be visually identical to original hgt when I compare them in QGIS.
First two methods I mentioned make TIFF files which has very light visual shift (~0.5 pixel).


texts:
I'm writing a script that calls a function I wrote in another script. I can't access the env path python is installed on so I must save my scripts in another location. To get one script to call another, I have to use sys.path.append(r"pathToPythonFiles") and then import the script as a module and then call the function.

Main script:

import sys
import arcpy
import os

sys.path.append(r"N:\Common\Script\Drivers")

import CO_181_ExStMetal as sMetal

sMetal.exportMetal()

script with callable function:

def exportMetal():

    # Create program variables

    # workspace
    mxd = arcpy.mapping.MapDocument(r"C:\RS_Data\Workspace\Batch_Carts_CO_181\Batch_Carts_CO_181.mxd")

    # dataframe
    df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]

    # field for SearchCursor
    numField = 'Route'

    # feature layer for cursor
    lyr = arcpy.mapping.Layer(r"C:\RS_Data\Workspace\Batch_Carts_CO_181\gisdb\layers.gdb\Service_Location")

    # create empty set to hold unique route values
    rtVals = set()

    # Create searchCursor
    cursor = arcpy.SearchCursor(lyr)
    for row in cursor:
        rtVals.add(row.getValue(numField))

    scenario = max(rtVals)

    path = r"S:\Shared\RouteSmart\Automated_BIC_and_Carts\CO_181\Standish_Yard\Metal"

    outNameCarts = os.path.join(path, ("Metal_Report.pdf"))
    print outNameCarts

    rlfCarts = r"N:\Common\Report Templates\CO_181_Metal.rlf"

    reportCarts = arcpy.mapping.ExportReport(lyr,rlfCarts,outNameCarts, "USE_RLF")

if __name__ == '__main__':

    exportMetal()

Here is the error:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec cmd in globals, locals
  File "N:\Common\Script\Drivers\sys_test.py", line 1, in <module>
    import sys
  File "N:\Common\Script\Drivers\CO_181_ExStMetal.py", line 10, in exportMetal
    import arcpy
NameError: global name 'arcpy' is not defined

I've tried import arcpy at the top of script 2 before defining exportMetal() as a global and then locally within the function as well. I also tried after if name=='main':

I'm concerned that once I set the system path in the script the program can't find arcpy which is installed in the python env path.


texts: Mapbox GL JS default state of toggle layers

texts:
How do I change the default state (off or on) of a toggle menu in Mapbox GL JS? The code that I'm using came from another question on this site Mapbox GLJS - group layers

<!DOCTYPE html>

    <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.20.0/mapbox-gl.css' rel='stylesheet' />


texts: Google earth incorrectly rendering coordinates?


texts: How to project ASCII grid given in degrees?

texts:
I have an ASCII file (.asc format) with 1/12 degree grid-squares of the whole world (see http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/), which I want to project onto WGS_1984_UTM_Zone_43S in order to determine some areas of land cover within some polygons in this projection.

ncols         4320
nrows         2160
xllcorner     -180
yllcorner     -90
cellsize      0.0833333333333333
row 1

I'm using ArcGIS at the moment, but when I put the ASCII straight into ArcMap, it incorrectly interprets the units as metres, not degrees (so the world map is only 360m from side to side!). This I put down to having no defined spatial reference

How should I set the spatial reference in order to get the appropriate scale?


texts: Fully cover polygon with circles of fixed radius?

texts: I have some detailed polygons I want to cover by generating points buffered to 4km radius. How do I determine the best way (or even a decent) way of covering the polygons completely with minimal overlap and fewest circles? What would be a good place to get started?


texts: Converting coordinate obtained from google maps to be used with gdal_translate in java or cs2cs?

texts:
this is a follow on question from Geo-Referencing an image from google maps with gdal_translate

I am trying to get cs2cs to work but I am getting the following error.

My command is as follows:

cs2cs +init=epsg:4326 +to +init=epsg:3857 -f "%%.2f" <wgs.txt >merc.txt

And the output when using PROJ_DEBUG=5 is

pj_open_lib(epsg): call fopen(C:\Users\monke\Desktop\Trial\cs2cs\Bin/epsg) - succeeded
pj_open_lib(epsg): call fopen(C:\Users\monke\Desktop\Trial\cs2cs\Bin/epsg) - succeeded
pj_open_lib(@null): call fopen(C:\Users\monke\Desktop\Trial\cs2cs\Bin/@null) - failed

I have been at this for hours and haven't made any progress.

Failing a cs2cs solution, a solution using java would be great too.


texts: GeoServer WMS OGC-filter fails with "I don't understand the tag: gml:Polygon"

texts:
I'm trying to construct an OGC filter to pass to a WMS GetMap request to GeoServer but I can't seem to do it. The goal is to show the parts of the layer that intersects with a supplied polygon.

When I run the query I don't get a result but get this in the server logs:

org.xml.sax.SAXException: Attempted to construct illegal filter - I don't understand the tag: gml:Polygon.  HINT: tags are case-sensitive!

Any ideas?

The query:

http://kar-toffel1.kartena.intra:8080/geoserver/kartena/wms?LAYERS=kartena%3Aosakert-luftkabel&STYLES=&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&SRS=EPSG%3A2400&FILTER=%3CFilter%20xmlns%3Agml%3D%22http%3A%2F%2Fwww.opengis.net%2Fgml%22%3E%3CIntersects%3E%3CPropertyName%3Ethe_geom%3C%2FPropertyName%3E%3Cgml%3APolygon%20srsName%3D%22http%3A%2F%2Fwww.opengis.net%2Fgml%2Fsrs%2Fepsg.xml%232400%22%3E%3Cgml%3AouterBoundaryIs%3E%3Cgml%3ALinearRing%3E%3Cgml%3Acoordinates%3E1271018%2C6404272%201271210%2C6404272%201271210%2C6404209%201271018%2C6404272%3C%2Fgml%3Acoordinates%3E%3C%2Fgml%3ALinearRing%3E%3C%2Fgml%3AouterBoundaryIs%3E%3C%2Fgml%3APolygon%3E%3C%2FIntersects%3E%3C%2FFilter%3E&CRS=EPSG%3A2400&BBOX=1270899.0723165,6404076.7312769,1271329.8765605,6404348.0392931&WIDTH=524&HEIGHT=330

The filter unescaped:

<Filter xmlns:gml="http://www.opengis.net/gml">
    <Intersection>
        <PropertyName>the_geom</PropertyName>
        <gml:Polygon srsName="http://www.opengis.net/gml/srs/epsg.xml#2400">
            <gml:outerBoundaryIs>
                <gml:LinearRing>
                    <gml:coordinates>
                        1271018,6404272
                        1271210,6404272
                        1271210,6404209
                        1271018,6404272
                    </gml:coordinates>
                 </gml:LinearRing>
            </gml:outerBoundaryIs>
        </gml:Polygon>
    </Intersection>
 </Filter>


texts: Polygon to polyline conversion in ArcMap: multiple entries in attribute table for the same line segment, why?

texts:
I have converted a polygon shapefile to polyline. I am having hard time understanding the polyline's attribute table. Below is part of the table:
objectid    left_fid    right_fid   shape_leng
2              44           6           .0117896
3              44           44          .0117896
4              6            6           .0117896
5              6            44          .0117896
6             -1            44          .0667701
7             -1            6           .0667701

I am looking at the border segments: the variable shape_leng indicates unique segments of the border. For example, objects 2-5 all indicate a single segment, with polygons 44 and 6 across the border. Objects 6 and 7 indicate another segment of the border, again with polygons 44 and 6 across the border.
What I don't understand is why the attribute table includes several rows for the same line segment? And how I should read left_fid and right_fid?


texts: Move shapefile upper and lower corner points

texts:
I have a shapefile that has been generated by an internal tool and had missing shx and dbf files but I was able to generate it using the post in the below link:

http://geospatialpython.com/2011/11/generating-shapefile-shx-files.html

However, when trying to perform spatial join with a set of points it returns null.

Upon checking out and comparing with normal shp files I found that upper and lower left and right points are very large

Example upper left is

-327802.687478,+3530832.250432

Is there a way to change those boundaries so I can perform spatial join or reproject the points I wish to do spatial join

I am using gvSIG and QGIS.


texts: Geo-reference a raster having separate latitude and longitude files in MODIS

texts:
R generated file:

Python generated file:

I want to georeference using external lat lng files in a MODIS data and then project it using R to WGS 84 coordinate system like how it happens using VRT file in Python. The code I could get so far shows wrong result in R when compared to result in Python.

    library(gdalUtils)
    library(rgdal)
    fn<-'MOD11_L2.A2017080.0005.006.2017082010345.hdf'
    sds <- get_subdatasets('MOD11_L2.A2017080.0005.006.2017082010345.hdf')
    filename <- 'Path/file.tif'
    gdalinfo <- GDALinfo(fn, returnScaleOffset = FALSE)
    gdal_translate(sds[1], dst_dataset = filename)
    r <- raster(filename) #how to georefernce it using the lat long files in subdtaset of the hdf file
    wgs1984 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
    projection(r) <- wgs1984
    extent(r)=c(xmn=130, xmx=170, ymn=-50, ymx=-20) # dont want to specify it since there are many files
    writeRaster(r, file="New_geo.tif", format='GTiff', overwrite=TRUE)
    plot(r)
    summary(r)

#summary of file generated using R
Min.    270.16
1st Qu. 304.76
Median  308.64
3rd Qu. 310.92
Max.    321.80
NA's      0.00
class       : RasterLayer
dimensions  : 2030, 1354, 2748620  (nrow, ncol, ncell)
resolution  : 0.0295421, 0.01477833  (x, y)
extent      : 130, 170, -50, -20  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : Path\file.tif
names       : file
values      : 0, 65535  (min, max)

#The files generated using python and warping shows this result. File generated in python gives this result
        X1_LST
Min.      0.00
1st Qu.   0.00
Median    0.00
3rd Qu.   0.00
Max.    321.58
NA's      0.00
class       : RasterLayer
dimensions  : 1747, 2485, 4341295  (nrow, ncol, ncell)
resolution  : 0.01217369, 0.01217369  (x, y)
extent      : 133.8272, 164.0788, -43.75374, -22.48629  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : Path\1_LST.tif
names       : X1_LST


texts:
I'm unable to complete processing in the Processing Modeler due to the following GDAL Error:

2019-01-21T10:47:12     CRITICAL    Could not create raster output E:\QGIS\QGIS_Projects\GDALTEST: GDAL provider Cannot load GDAL driver
2019-01-21T10:47:12     CRITICAL    Error encountered while running Reclassify by table

Here are the settings:

Using QGIS version 3.4.3


texts: Intersect five different points of same layer

texts:
In ArcGis Desktop I have the following problem: I have a layer of sites, let them be hotels, for instance, and I want to find each place, where there are less than five hotels within a radius of 30 km.

I thought this might somehow work with iterating over the intersect function in the Modalbuilder. However, I have never used the iteration tool and am not quite sure, whether I'm doing this right. Additionally, how do I make sure, that it doesn't intersect with the same point twice, when I iterate over it?

Is there a different approach that would be better for this?


texts:
New to QGIS and really all things GIS.

I'm looking at specific types of attacks in a region. Surrounding these attacks are other points like schools, clinics, political buildings, etc. I'm trying to identify what is in each of the buffers (a radius around the attack location).

I'm able to do it for each variable but would rather do it as a whole for everything, i.e. across the 23 attacks, what important sites are within a 10km radius of each attack?


texts: Geotagging of batched photos with large size

texts:
I am trying to geotag some pictures in ArcGIS 10.5. I have pictures contain coordinates and want to place pictures on the final layout of the map. I have tried two methods but did not get the desired results.

I used Geosetter tool to convert pictures in KML format, after that I imported the KML file in arcmap and I got below result. But here picture size is too small to get any info from it. I just need the same output with the large size of pictures.

Secondly I used Geotagged photos to point tool of ArcMap. In this case I can't display pictures in layout view.

How can I display the pictures on final exported map as shown in first figure but with larger size?


texts: Using Python to count records and show difference in ArcGIS Desktop?

texts: I'm trying to find an example of a python script for ArcGIS 10.1 that will get count of a feature class (FC) in my ArcSDE and get count of a FC in an incoming file geodatabase update and then write the difference to a log file. We get monthly updates from a vendor. For instance, wells FC, and want to be able to compare record counts and log diff (a-b=C). Right now I have a script that will record start time, get count, truncate, append, rebuild index.


texts:
So I have been doing a lot of research on this topic but I am new to working with online mapping. I have an access database and I would like to portray information in a google map. I have learned how to embed a google map into an access database but I am looking for the reverse. I want to take the data from the access database export it out, send it to a website, and have a google map on the website showing the information.

Any suggestions on how to perform this action?


texts: ogr2ogr Shapefiles to Postgres SRS problem

texts:
I am trying to import to my Postgresql (9.1+129ubuntu1) + PostGIS (1.5.3-2) database a Shapefile with polygons with the ogr2ogr tool from Gdal-Bin package (1.7.3-6ubuntu3).

I know the origin SRS of the Shapefile (EPSG:23031), and I know the destination SRS used in the database so I can represent all the diferent objects together (EPSG:25831), so I use this command:

ogr2ogr -f "PostgreSQL" "PG:dbname=DBNAME host=127.0.0.1 user=DBUSER password=DBPASS" "/PATH_TO_SHAPEFILES/SHAPEFILE.shp" -s_srs EPSG:23031 -t_srs EPSG:25831 -lco PRECISION=NO -lco GEOM_TYPE=geometry -lco GEOMETRY_NAME=the_geom -lco EXTRACT_SCHEMA_FROM_LAYER_NAME=NO -nlt GEOMETRY -nln DBTABLENAME

And of course in the path "/PATH_TO_SHAPEFILES/" I have the files SHAPEFILE.shp, SHAPEFILE.shx, SHAPEFILE.dbf and SHAPEFILE.prj

But after the import, I try to represent that polygons and they seem moved to de East, so I guess is a projection transform problem or definition.

Here is an snapshot of the representation in the map by OpenLayers. The correct situation of the polygons should be the grey lines, and the imported polygons are the blue polygons with black lines:

EDIT: As @tilt recommended me I tried to do the same process but converting from shapefile to shapefile, with this command:

ogr2ogr -f "ESRI Shapefile" pruebashp2shp.shp "/PATH_TO_SHAPEFILES/SHAPEFILE.shp" -s_srs EPSG:23031 -t_srs EPSG:25831

I opened the result by QuantumGIS together with the original layer, and it is also moved to the East, as the following image shows:


texts:
I am creating a map feature editor that supports versioning & rollback using Geoserver as the WFS.

From what I've read online, it seems Geoserver supports feature versioning. My question is how do I implement versioning on Geoserver? Do I need to make any sort of modification on the data store for this to work or will the versioning be handled automatically by Geoserver?

Would be great if there was a tutorial on this.

Thanks.


texts:
I've imported some of line shapefiles drawn in AutoCAD, so now I have one shapefile with multiple features (one feature is one line drawn in AutoCAD).
How can I merge all of these features into one single feature?


texts: Converting projected geoTiff to WGS84 with GDAL and Python

texts:
Apologies if the following question is somewhat stupid, but I am only VERY new to this whole GIS thing.

I am trying to convert some projected geoTiff images to WGS84 using gdal in python. I have found a post that outlines the process to transform points within projected GeoTiffs using something similar to the following:

from osgeo import osr, gdal

# get the existing coordinate system
ds = gdal.Open('path/to/file')
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system
wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs)

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5]

#get the coordinates in lat long
latlong = transform.TransformPoint(x,y)

My question is, if I want to convert these points and create a new WGS84 GeoTiff file, is this the best way to go about it? Is there a function that exists that will do such as task in 1 step?

Thanks!


texts: Transferring shapefile to Garmin Monterra with labels/attributes

texts: I would like to transfer a shapefile to my Garmin Monterra, but preserve at least one field of attributes.  So far, I have tried converting to kml, but when I load into basecamp, there are no attribute data and the waypoints are all given generic numbers.  Is there a way to export shapefiles from ArcMap so that the loaded point data in basecamp can be labeled according to an original field of attributes?


texts:
I am trying to build an arcgis nearby places android quite similar to this:
https://developers.arcgis.com/example-apps/nearby-android/. I have make a map portal and import geojson file which placed all my stations location.

But i can't get the app to show my location.

package com.utbfyp.webmapapp;

import android.support.v7.app.AppCompatActivity;
import android.os.Bundle;
import android.app.Activity;
import android.location.Location;
import android.location.LocationListener;
import android.view.View;

import com.esri.arcgisruntime.mapping.ArcGISMap;
import com.esri.arcgisruntime.mapping.view.LocationDisplay;
import com.esri.arcgisruntime.mapping.view.MapView;
import com.esri.arcgisruntime.portal.Portal;
import com.esri.arcgisruntime.portal.PortalItem;

public class MainActivity extends AppCompatActivity {
    private MapView mv;

    @Override
    protected void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        setContentView(R.layout.activity_main);
        mv = (MapView) findViewById(R.id.map1);

        LocationDisplay myLocation = mv.getLocationDisplay();
        myLocation.startAsync();

        //use ArcGIS Online as the Portal
        Portal arcgisOnline = new Portal("http://www.arcgis.com");
//create a PortalItem with ArcGIS Online WebMap ID
        PortalItem portalItem = new PortalItem(arcgisOnline, "xxxxx");
//create an ArcGISMap object with the PortalItem
        ArcGISMap myMap = new ArcGISMap(portalItem);
//display the webmap to MapView
        mv.setMap(myMap);

            }

        }


texts: Turning sublayers on/off in QGIS2Web Leaflet export

texts: For some simple maps, I need to turn sublayer on and off. Like in the example (a) and (b) with Leaflet/Qgis2Web.


texts: Difference between Geographic and Projected coordinate systems?

texts: What is the difference between a geographic coordinate system and a projected coordinate system?


texts:
What are the valid extents of the British National Grid?
Are all tiles AA through ZZ valid, or is it just those which start with H,J,N,O,S,T that should be used?


texts: Error on Reconcile - Inconsistent Network State

texts:
I am encountering an error that I am not sure how to address when I reconcile an edit session.  I'm using ArcGIS 10.2.1, ArcFM 10.2.1a, connected to a SQL 2012 database via a direct connection using operating system authentication on a 64 bit Windows 7 machine.

I am receiving the error while attempting to reconcile a session containing a large number of edits removing loops from an electrical network.

Here is the error I receive:

I think I have a good understanding of what the error is describing - I just can't figure out how to address it.

My first thought was to add the Network Junctions to the map and look for junction 85752 and 90717 but when I add the junctions it comes in from DEFAULT, not the version that I'm trying to reconcile and those junctions don't exist in default (or at least I can't find them in the table).  So that seems to be a dead end.

I've also tried running the "Repair Connectivity" and "Verify Connectivity" tools from the Geometric Network Editing toolbar (both were run on the entire network - NOT the current extent) in hopes that those would deal with the inconsistency and then it would reconcile.  Both tools ran without incident and neither had any impact on the result.

My fundamental question is: Is there anything I can do about this or do I just need to delete the version and start over?

--EDIT--

I found an ESRI help article that addresses this issue (I've been through the ESRI help article on this error (http://support.esri.com/es/knowledgebase/techarticles/detail/21187) it seems like my only option if the tools don't run right is to delete the version and start over.  Which isn't great news but I thought I'd post it here in case someone was curious.


texts:
I'm using QField to collect geologic data in the Outback.

I'm creating point data, with various fields for rock type etc.

Occasionally I take samples of the rocks, and want a field that auto-increments the sample number. The problem I have is that most entries do not require a sample, so a simple '$rownumber' query is not suitable.

What I'm picturing is a field for "sample_taken?"   'Y'/'N', and if 'Y' then an auto-incremented number is inserted continuing from the last sample number I used. For example:

Thanks in advance!


texts: geoserver wrong tile generation on leaflet

texts:
I am creating world map with bathymetry using Geoserver 2.5.1 and leaflet.
It works perfactly upto certain zoom level. but after that i get few tiles drawn in wrong color as shows in images. Does anyone have any idea why geoserver behaving like it or any direction to fix it?

As i zoom in more on the map tiles are missing colors as shown below.


texts: Convert Natural Earth shape to GeoJSON - (Multi)LineString instead of (Multi)Polygon?

texts:
I need shapes of admin regions as GeoJSON (Multi)Polygons. After looking for an appropriate source for way too long, I found this NE Admin 1 dataset.

Using ogr2ogr to convert the data to GeoJSON, I found features with geometry-type (Multi)LineString. Not quite what I wanted. This guy on GitHub used a different dataset from Natural Earth, but he got (Multi)Polygons. Are there some special options for ogr2ogr, some possible conversions, or is it the dataset?

I need to determine whether a point is in a province/region and with the tools I'm using, that's easiest with polygons. If it happens to be a problem with the dataset, does anyone know an alternative?

That's the command I used for conversion:

ogr2ogr -f GeoJSON -t_srs crs:84 output.geojson ne_10m_admin_1_states_provinces_lines.shp`

PS: Most examples I've found specify -t_srs crs:84. I do need WGS84, as I'm also using data from geonames.org. But I'm not sure if I even need the transformation.


texts: esri-leaflet edit multiple featureLayers

texts:
I'm trying to use the esri-leaflet example on editing (https://esri.github.io/esri-leaflet/examples/editing.html) to edit two different featureLayers.

First when I add a point with leaflet draw, I want it to automatically create a buffer around the point. I can get it to add the circle in 'draw:created' by pulling the points (lat & lng) from the marker and giving it a radius.

 map.on('draw:created', function (e) {
    // add the feature as GeoJSON (feature will be converted to ArcGIS JSON internally)
    collectionBoxes.addFeature(e.layer.toGeoJSON());
   //  createBuffer using the lat lng from marker added by leaflet draw
    var circle = L.circle([e.layer._latlng.lat, e.layer._latlng.lng], 304.8).addTo(drawnItems);
   // try to add buffer to arcgis service
    buffers.addFeature(circle.toGeoJSON());
    disableEditing = false;
});

This line adds a circle on the map.

  var circle = L.circle([e.layer._latlng.lat, e.layer._latlng.lng],
  304.8).addTo(drawnItems);

This line does not give any errors.

  'buffers.addFeatures(circle.toGeoJSON);'

but does not edit the arcgis service for the buffer. The point layer (collectionBoxes) does get edited. Any ideas on this or how to edit the example to edit two layers?


texts: Unable to delete GeoTIFF after adding it to GeoServer

texts:
After adding a file to GeoServer via REST I am unable to delete it from the directory on my computer. I always get the same error saying that the file is in use by another process.

This only happens when apply a style to it that will show the layer as a contour map.

I have tried removing the layer from GeoServer but the GeoTIFF will remain in use anyway.


texts: Drawing equal paths from geographic points?

texts:
I have 10'000 geographic points (POS), 10 salesmen, and one starting point. To avoid discrimination between salesmen, everyone should travel the same distance and the same number of points of sale.

My goal is to:

assign a path for each salesman starting from the starting point with the same number of geographic points (POS).
Each path should have the same distance.

My solution is sort an array of GPS coordinate points by distance and then assign each point to a path until the end of the loop, but doubt this will work..

What is the easiest algorithm to achieve this?


texts: How to convert shapefile to a set of Latitude and Longitude coordinates in QGIS?

texts:
I am doing research in water and climate. I have downloaded a map for World Water Bodies that describes all the large water bodies that are drawn as polygons, and another map for small water bodies (World Linear Water)that are drawn as lines.

I am not a GIS expert, and I have only QGIS software, so I followed this post to load the LPK file to QGIS then I tried this post and this post to convert this layer to a point-set of Longitude and Latitude coordinates. Unfortunately, I stuck in adding the calculated fields $x and $y to the layer, as it lasts for hours without real progress but "Not responding" while the qgis-bin.exe process in the Windows Task Manger shows about 2-10GB memory and 14% CPU of a core-i7/16G-RAM PC.

Layer 1 is 528MB size, and the Nodes file is 3GB with 38 Milion features!

Is there a method to speed up the conversion?

Am I doing the correct steps?

I just want a Long/Lat list (CSV) of the points surrounding the polygons of continents, seas, lakes, and big rivers, to use them in my research.


texts: Automatic matching of two street networks

texts:
I'm faced with the problem of matching two street networks. I have these two differently sourced networks. The links in each network can be segmented, i.e. there is not necessarily a single whole link between two intersections, instead the link can be broken into an arbitrary number of segments with different attributes.

I should transfer the attributes from one network to the features of the other network and I figured for doing that I will have to match corresponding streets/street segments.

My current idea is to densify the links in either network and then do an analysis of the cumulative distances between vertices in the two networks? If the nodes of two segments overall have small distances from each other, the two links probably correspond? However, I can see problems if two corresponding road segments are expressed in links of wildly different length in the two networks. That notwithstanding, I'm unsure as to how to best go about implementing a method that would analyse the two networks and maybe also offer a way of visually exploring the quality of the matching.

On the plus side I'm familiar (to different degrees) with many platforms (ArcGIS, QGIS, SAGA, R, Python, ...) and can install any promising tool. I appreciate any pointers you may be able to give.


texts: I'm working on building a GIS for Electric power distribution system. Can someone help me on representing an electric switch( With On/OFF state which will connect /disconnect downstream lines from the network) in QGIS. I want to detect the number of customers who will get affected due to turning off an upstream switch in a power line.


texts: What programming languages are used by different GIS platforms?

texts:
Suppose a programmer has no knowledge of GIS, but wants to leverage his/her programming skills to get into GIS.  Please list different GIS platforms, packages, etc. that use a particular language.

(I was thinking, as suggested, about doing a simple test of the code formatting capability of StackExchange.  I figure this might be done in the context of a useful answer.)

When you answer, if the language is one of those listed below, please post a code snippet and test the formatting to make sure it works.  Here's an excerpt from Jeff Atwood's answer on meta-stackoverflow.

The spec is:

<!-- language: lang-or-tag-here -->

    code goes here

<!-- language: lang-or-tag-here -->

    code goes here

You may use either a tag or a prettify language code to specify, though prettify language codes are always guaranteed to work regardless of what language the tag happens to be set to.

Use <!-- language-all: lang-or-tag --> to use the given highlighting for all the following code blocks.

The complete list of supported prettify code languages is:

lang-none
lang-default
lang-bsh
lang-c
lang-cpp
lang-cs
lang-csh
lang-css
lang-hs
lang-html
lang-java
lang-js
lang-lisp
lang-lua
lang-ml
lang-perl
lang-php
lang-py
lang-proto
lang-rb
lang-scala
lang-sh
lang-sql
lang-vb
lang-xml


texts: How to join Shapefiles without common attributes using ArcMap

texts: I have a Dissemination Areas boundaries shapefile with has an attribute table full of data needed for lowest quartile analysis and I have a neighborhood shapefile which has lowest income data which also needs to be considered but am not sure how I can combine these shapefiles so I can consider both the lowest income neighborhood data as well as the dissemination area data in my analysis of worst overall dissemination areas. I tried spatial joining but they don't share a similar value.


texts: How can my PyQGIS script copy specified attributes from one layer to another?

texts:
I am trying to create a script that will take selected features from one layer and search for the specified attributes related to a column name. in my case 'id'. What I am having trouble doing right now is writing the script to tell qgis to take the selected 'id''s and write them to a new attribute table under a specified field.
Below is my current script.

from qgis.gui import QgsMapTool
print 'imported'
iface.showAttributeTable(iface.activeLayer())
print 'Attribute table shown'

# define variables
layer = iface.activeLayer()
feature = layer.getFeatures()
selection = layer.selectedFeatures()
# prints the number of attributes the user selected.
print len(selection)

# for loop to identify trail_stuctures from attribute table
for feature in selection:
    print feature['id']
output_file = open('c:/Users/Pulse3/Desktop/plugin.txt', 'w')

# Writing selection id's to a text file stored on the desktop
for selection in layer.selectedFeatures():
  line = '%s\n' % (selection['id'])
  unicode_line = line.encode('utf-8')
  output_file.write(unicode_line)
output_file.close()

The above script is only one of the two scripts. Right now the user has to click on the layer through the layer panel. manually go into the attribute table, and use the select feature tool to select a series of lines. ( which are my spans.) Only the selected lines are used in this script.

Ultimatly id like the script to be able to look at seleced features in one layer and Write the 'id' into a standalone table. then switch to another layer and look at the user selected features of another layer and do the same thing.

Summary:

Two layers:

Span

Service_point

Both layers will have a selection done based on a user selected features with the mouse.

Both layers will also need to be able to take the returned attribute values list and write it into two different field within a seperate stand alone table.

Essentially I am designing fiber internet drop maps and we may have service points which are houses. that each have a unique ID number and then the spans which are lines off the main backbone fiber line that also have unique id's. both of these features are in their own layers.

I am trying to write a script that allows me to extract just the Id from the service point and then the user selected spans and fills in new geometry attribute table with the same service point id for each of the spans that feed the service point and then the span id for each of the span ids that service the service points.

Is there a way to do this in QGIS?


texts: Function (sample code) to extract raster value per polygon in R

texts:
I have a raster image (Sentinel-2 band 4) and a shapefile covering the same area than the image (same projection and extension). The shapefile has different polygons which are agricultural fields.

I want to create a new file where, for each polygon, the mean value of the pixels that are inside it is calculated.

Could you please help me to find an appropriate R code for that?


texts: Extract elevation by longitude and latitude

texts:
I have a large file of coordinates on the map of New Zealand, by longitude latitude. I want to find the approximate elevation at each point.
Starting with just one area, I have found a raster with the data I need. I have loaded it into QGIS, and it looks good to me. It has this information with it:
Name    NZDEM_SoS_v1-0_27_Dunedin_gf
Path    C:\...\elevation\kx-27-dunedin-15m-dem-nzsosdem-v10-GTiff\NZDEM_SoS_v1-0_27_Dunedin_gf.tif
CRS EPSG:2193 - NZGD2000 / New Zealand Transverse Mercator 2000 - Projected
Extent  1372000.0000000000000000,4866000.0000000000000000 : 1492000.0000000000000000,5046000.0000000000000000
Unit    meters
Width   8000
Height  12000
Data type   Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF

So far, so good. Now in R, I have installed package raster, and I can do:
fname = "../elevation/kx-27-dunedin-15m-dem-nzsosdem-v10-GTiff/NZDEM_SoS_v1-0_27_Dunedin_gf.tif"
elev.r <- raster(fname)

This gives me this warning, but perhaps it isn't a problem:

Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
but +towgs84= values preserved

Then I can do
extract(elev.r,1000,1000)

and that returns a value 1125.455  which is presumably the elevation.
What do I need to do to convert my longitude, latitude to the x,y which the extract function will understand?
I downloaded the raster here: https://koordinates.com/my/downloads/2000967/download/?dl
long = 170.605375
lat =  -45.859668

xy <- cbind(lat,long)
colnames(xy) <- c('x', 'y')
xy <- as.data.frame(xy)

coordinates(xy) <- ~ x + y # telling R these are spatial points
crs(xy) <- crs(elev.r) # set the same crs as in your_raster
crs(xy)

extract(elev.r, xy)

returns NAs whatever appears to be in xy


texts: GetLegendGraphic return a symbol instead of colormap ramp

texts:
I'm trying to put inside an OpenLayers project the legend for an NDVI map. Previously I've created with QGIS 3.10 the SLD as indicated on the GeoServer's official documentation.
When I try to use GetLegendGraphic to see the legend's preview I see only the symbol of raster and not the legend. If I set the style as default for my raster layer I can use the link from GetLegendGraphic for see correcty my legend. But If I set "Associated" instead of "Dafault" I see again the raster's symbol.
Why I'm obligated to use "Default" for styling a raster?
EDIT
Here the WMS stylized(A) and here the WMS not stylized(B). The style name is NDVI.
https://gis.massimilianomoraca.it/geoserver/MassimilianoMoraca/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetLegendGraphic&FORMAT=image%2Fpng&LAYER=NDVI_Campania_20150807&SCALE=136494.42037699942
gives

it is possible to see the legend and is not possible to do the same with
https://gis.massimilianomoraca.it/geoserver/MassimilianoMoraca/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetLegendGraphic&FORMAT=image%2Fpng&LAYER=NDVI_Campania_20160712&SCALE=136494.42037699942
.
The layer A has NDVI as default style and for B I've used "Associated".


texts: Issues with GeoServer GetFeatureInfo request parameters

texts:
I managed to make GetFeatureInfo request but it has some issues still.
This is code:
function request(e) {
  var BBOX = that.mymap.getBounds().toBBoxString();
  var WIDTH = 101;
  var HEIGHT = 101;
  // var X = that.mymap.latLngToContainerPoint(
  //   e.latlng,
  //   that.mymap.getZoom()
  // ).x;
  // var Y = that.mymap.latLngToContainerPoint(
  //   e.latlng,
  //   that.mymap.getZoom()
  // ).y;
  var X = 50;
  var Y = 50;
  let URL = `/geoserver/cite/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=true&QUERY_LAYERS=cite:tbilisi3&STYLES=&LAYERS=cite:tbilisi3&exceptions=application/vnd.ogc.se_inimage&INFO_FORMAT=text/html&FEATURE_COUNT=50&X=${X}&Y=${Y}&SRS=EPSG:4326&WIDTH=${WIDTH}&HEIGHT=${HEIGHT}&BBOX=${BBOX}`;

  console.log(URL);
  $.ajax({
    url: URL,
    datatype: "html",
    type: "GET",
    success: function (data) {
      console.log(data);
      // L.popup().setLatLng(e.latlng).setContent(data).openOn(that.mymap);
    },
  });
}

I have two problems:

You can see the bbox contains all visible map, hence it makes request on larger region, how do I make this region smaller? e.g. when user clicks I only want info from a smaller region around the click, not whole visible map.

This request works as it is now just there are many results due to bbox. But, if I uncomment X and Y and use them (instead of hardcoded values), the response doesn't arrive anymore and I think I also get error about invalid point. What is wrong with  my (commented) X and Y parameters?


texts: Changing colors in SVG file polygon using QGIS?

texts:
I try to change colors and lines width in this polygon (this is SVG file):

I learn from here
SVG symbol colors differences in QGIS 2.12?
Using SVG file to fill polygon using QGIS- problem
to add this line to the SVG file in order to change colors:
fill="param(fill) #FFF" stroke="param(outline) #000" stroke-width="param(outline-width) 1"

but the option to change the fill is blocked

How can i change the fill color?


texts: Drop-down list from Existing Column to Select 1 Record to use in Model?

texts:
It has been a while since I've used ModelBuilder or Python.

I have a point feature class
I want to be able to use a drop down list to select one of those points
Drop down will be based on a Name field in the point feature
That point will then be fed into the rest of my model

Is there a way to do this with ModelBuilder/tools, or just as a simple python script insert into ModelBuilder?


texts: How to delete a bracket in the attribute table of QGIS?

texts:
I would like to delete a bracket ([]) in one column of my attribute table. I have different numbers in the column, written like this "[901], [555]..." How can I delete the bracket, withouth having to edit each field?

I am using QGIS version 2.12.3-Lyon.


texts: How to get travel time and link used from osm2po http service

texts:
I used osm2po to get shortest paths but I would like to get travel time and osm links ids used in the path. Is it possible?

Here is the osm2po service url format I use (this is an example):

http://localhost:8888/Osm2poService?cmd=fr&source=26915&target=106044&findShortestPath=false&ignoreRestrictions=false&ignoreOneWays=false&routerClassId=0&heuristicFactor=0.0&maxCost=0.0&debugExpansion=false&key=Key&value=Value&tsp=

Thanks!

Leo

Comment:
Well, I just found I can add format=geojson to the query string and I can then get the line strings associated with the path. Perfect. I then just need to match the id in the db with the travel speed.


texts: Overpass output city name instead of coordinates?

texts:
I did a query to look for certain things. However it only outputs the lat and long coordinates, there is no city data available. I'd like to have the city it belongs to / closest city though.

I know there are websites where I can convert a single coordinate to a city.
But since my initial output consist of hundreds of elements, is there a way to alter the query to convert it to the closest city in the first place?


texts: Image orthorectification with XYZ and roll, pitch and yaw

texts:
There are a lot of raw images taken by UAV and the only parameters that I have are XYZ and roll, pitch, yaw in excel spreadsheet.
I've been trying to orthorectify these images with ArcGIS, Er mapper and even Python. They all require other paramteres such as control points, camera files, ... Which unfortunately I don't have them.
Is there a way to orthorectify these images using only XYZ and roll, pitch, yaw?


texts: Edit Field Value of Feature Class Using Python in ArcGIS 10.2

texts:
I have a feature class whose attribute table looks like this:

Latitude | Longitude  | Area

For the area column, all rows are null.  I would like to put in a specific value for all rows using Python.

I know how to do this on the front-end in the ArcMaps application ( right click > Open Attribute Table > right click on column > field calculator > input value).

How can I accomplish this programmatically using Python?


texts: How to get WorkspaceCategory for a given path

texts:
I need to determine the WorkspaceCategory string (IGxDialog.FinalLocation.WorkspaceCategory) for a given path. I need this to create appropriate WorkspaceFactory and finally be able to create multiple feature classes / shapefiles in a given location.

Right now I'm using IGxDialog with GxFilterContainers to select final location, but when I select a file geodatabase, the pGxDialog.FinalLocation.WorkspaceCategory returns "Folder" (the parent folder of the selected container), not the category of the selected container.

    Dim pGxDialog As IGxDialog = New GxDialog
    Dim pFilterCol As IGxObjectFilterCollection
    Dim pGxFilter As IGxObjectFilter
    Dim pEnumGx As IEnumGxObject
    pGxDialog.AllowMultiSelect = False
    pFilterCol = pGxDialog
    pGxFilter = New ESRI.ArcGIS.Catalog.GxFilterContainers
    pFilterCol.AddFilter(pGxFilter, True)

    If (pGxDialog.DoModalOpen(0, pEnumGx) = True) Then
        sOutputPath = pGxDialog.FinalLocation.FullName
        sOutputName = pEnumGx.Next.Name
        workspaceCat = pGxDialog.FinalLocation.Category
    End If

So, the questions are:

how can I get the WorkspaceCategory string for a given path?
are there any other ways of creating appropriate WorkspaceFactory for a given location (e.g. if my output location is a folder, I would need to create ShapefileWorkspaceFactory, if the output location is a file geodatabase, I would need a FileGDBWorkspaceFactory etc.).


texts: How to cut pixel size rectangle from GeoTIFF image with GDAL Python

texts:
I have GeoTIFF image with EPSG "4326" and longitude, latitude of a point in the image.

I need to cut a sub-image from the original image around the point (lng,lat) with the size of 128x128 pixels.


texts: Deleting duplicate records in the attribute table and keep the minimum value in PyQGIS?

texts:
I have two fields: ID_PARCEL and NL. The ID_PARCEL is the field that contains duplicate records. I would like to remove the duplicate record by keeping the minimum NL field value. For example: the value 22 in the ID_PARCEL field is duplicated, so the row of the attributes table that should be removed is the one with the highest value in the NL field, which in this case is the value 627.

I searched on the site, but what I found was deleting duplicate records without the condition to keep the minimum value of another field. Can someone help me?


texts: Remove a QgsRectangle from a RectangleMapTool in PyQGIS

texts:
In a QGIS plugin, I implemented a RectangleMapTool that I apply to my canvas. Thus, the user has the possibility to select a bounding box to define the spatial processing footprint.

Once he has selected his bounding box, I'd like to delete this QgsRectangle but I can't do it. If the user clicks anywhere after drawing the bounding box, the rectangle disappears but sometimes it remains.

You can run the following code in the Python console. If you run the script a second time, the bounding box will remain drawn:

from processing.gui.RectangleMapTool import RectangleMapTool

# Function
def bbox():
    r = tool.rectangle()
    line_bbox.setText(str(r.xMinimum()))
    line_bbox.show()

# MapTool
canvas = iface.mapCanvas()
tool = RectangleMapTool(canvas)
tool.rectangleCreated.connect(bbox)
canvas.setMapTool(tool)

# QLineEdit
line_bbox = QLineEdit()
line_bbox.hide()

I explored the documentation to be able to delete a QgsRectangle. When it is an entity of a layer I know how to do it but in my case, it is drawn on the canvas.

I also explored several links that ran iface.mapCanvas().scene().items(). I tried to remove the QgsRubberBand objects but my QGIS crashes or no changes are made.


texts: GFS: Global Forecast System 384-Hour Predicted Atmosphere Data, making daily precipitation and temperature forests

texts:
I am trying to understand and use the "GFS: Global Forecast System 384-Hour Predicted Atmosphere Data". As I got it, the dataset consist of forest for coming 384 hours (16 day) in hourly and 3 hourly scale. So first 120 hours in hourly scale then 3 hourly scale. Please, correct me if I am right and help me to figure out. Using this dataset I would like to make a code which will make a chart with forested daily rainfall, max temperature and min temperature. And I have tried to develop a code for daily precipitation, but does not work. Please, help me to improve it.

var now = ee.Date(new Date().getTime());
print(now);

var dataset = ee.ImageCollection('NOAA/GFS0P25')
                .select('total_precipitation_surface')
                .filterDate(now);
print(dataset);

var day1_pcp = dataset.filterDate(now,now.advance(23,'hour')).sum();
var day2_pcp = dataset.filterDate(now.advance(24,'hour'),now.advance(47,'hour')).sum();
var day3_pcp = dataset.filterDate(now.advance(48,'hour'),now.advance(71,'hour')).sum();
var day4_pcp = dataset.filterDate(now.advance(72,'hour'),now.advance(95,'hour')).sum();
var day5_pcp = dataset.filterDate(now.advance(96,'hour'),now.advance(120,'hour')).sum();

var days = dau1_pcp.merge(day2_pcp).merge(day3_pcp).merge(day4_pcp).merge(day5_pcp)
print(day1_pcp);

var chart = ui.Chart.image.series(summedCollection, wua_filtered, ee.Reducer.mean())
  .setChartType('ColumnChart');

print(chart)

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


texts:
I'm in need of listing all features in a linestring layer that have dangles (that is, one or both of its boundaries do not touch another feature). This being done in GDAL/OGR.

What I'm doing is iterating feature by feature, getting the bounding box of it, using that to do a spatial filter so I only select the features close to the one I'm checking, and then doing topology checks between it and the selected features. Something like this:

for i in range(n_feats):
    feat = layer.GetFeature(i)
    # get start and end points of feature
    # get bounding box of feature

    layer.SetSpatialFilterRect(bbox[0], bbox[3], bbox[2], bbox[1])
    n_feats_select = layer.GetFeatureCount()
    for dummy_i in range(n_feats_select):
        feat_select = layer.GetNextFeature()
        # check if feature intersects with previously gotten start and end points
        # if they don't, feature has a dangle, do stuff

    layer.SetSpatialFilter(None)

Now, this script works, which certainly matters, but it's slow. Specifically, I'm dealing with a layer with about 2 million features, and the spatial filter is hogging a lot of overhead. Specifically, each feature is taking about 10 seconds to process (the filter method alone taking over 9 seconds), thus 2 million features at this pace would take 231 days to process, give or take.

Is there any GDAL method that already does this, or a better algorithm for accomplishing this task?


texts: Adjusting same CRS to overlay raster satellite image with vector data

texts: I need to overlay my raster image of chollian satellite with shape file of South Korea. It worked in MODIS images but it does not work in chollian satellite images. Could you let me know how to solve this problems?


texts: Point style with height, width and angle

texts:
I'm searching for a special style for a point shape. In attribut table, I have, among others, 3 fields: height, width and angle.
I would like to use these 3 fields to represent my points as turned rectangles.

Somebody have an idea to do this ?


texts: GeoDatabase (GDB) import not working (QGIS3)

texts:
I have a *.gdb which I normaly could load in QGIS version (3.4 / 3.8) following that instructions (or symply drag'n'dop): Installing File Geodatabase (*.gdb) support in QGIS?

But it returns a error: Invalid Data Source: [...] is not a valid or recognized data source.

I also tryed ogr2ogr to translate it first in to shapefiles:

ogr2ogr -f "ESRI Shapefile" ./shp ./Schoeningstedt.gdb

With no result.
Whats wrong? Could anyone test my gdb? You can download it here: www.map-site.de/tmp/Schoeningstedt.gdb.zip


texts:
Another rasterio follow-up question, after questions on reading and writing rasters. I've successfully made and reprojected an output layer, but while trying to use the overview module to build pyramids, I found I don't understand something fundamental: the ctx parameter of the function.

def overview(ctx, input, build, ls, rebuild, resampling):

By a Google search, I found that ctx means "context" and that there's a module named context in click library used by rasterio. I got stuck here: further searching shows "context" only in other contexts. The only useful hints I found here: info.py and test_rio_info.py files in the rasterio repository. I'm not good enough to be able to reverse engineer the ctx format expected by the overview function from these libraries (not sure whether it is possible).

So I need a good manual explaining contexts as used in rasterio, either as an answer or as a link with some summary. How to produce/get them, what exactly do they mean (do they have any use besides passing data into functions?), what are their components (if any) etc.

EDIT: I've found the Context class in the click library in the core.py file on my disk (/usr/local/lib/python2.7/dist-packages/click/core.py on my Ubuntu machine). I didn't manage to explicitely initialize this class object (I didn't understand the params), but after reading Kersten's answer I think it is not needed while using Rasterio through Python API. The class documentation seems to agree:

  The context is a special internal object that holds state relevant
  for the script execution at every single level.  It's normally invisible
  to commands unless they opt-in to getting access to it.

  The context is useful as it can pass internal objects around and can
  control special execution features such as reading data from
  environment variables.


texts: Adding CRS definition/details to layout in print composer

texts:
I'm new to QGIS and I'm wondering if there's any way I can add a label to my map in composer with a similar format shown in the image below. I know that I can add the CRS and CRS definition using the layer property format but I would like a similar output to the one in the image.

The attached image is generated from Arc by the way.


texts:
Looking at the QGIS API Documentation:

bool QgsFeature::hasGeometry
    Returns true if the feature has an associated geometry.

Is invalid geometry different from a feature having no associated geometry? If so, is there a tool like .hasGeometry in PyQGIS that returns whether or not the geometry is valid?


texts: Incomplete type description in Geoserver WFS 1.1.0 and 2.0.0

texts:
I have a PostGIS database (2.1.3 on PostgreSQL 9.3.5) exposed as a WFS service with the geoserver (2.6-RC1).
This table has 4 columns: an id, a String, a Date and a Point Geometry.

When calling DescribeFeatureType with Version 1.0.0, the XSD describes a complexe type with the 4 columns.
But when I use 1.1.0 or 2.0.0 the XSD does not include the gml:PointPropertyType attribute refering my last colum of type Point Geometry.
I'm then not able to do use any geographic/geometric filter on that WFS service.

The only difference I've seen is that the GML version used in reference in the output XSD is not the same (GML 2.1.2 for WFS 1.0.0, and GML 3.x for WFS 1.1.0 and 2.0.0).

What did I miss to let this work in 2.0.0 or 1.1.0 at least?
PS: I need WFS > 1.0.0 to use time filters on my date column.

Edit:
When configuring WFS service inside Geoserver, it seems that the "Override GML Attributes" check-box for GML 2 what checked by default, while for GML 3 and 3.2 it wasn't. Once these checkboxes are checked the returned XSD also contains my fourth column of type gml:PointPropertyType.