texts: Why does this Geotools Java reprojection to WGS-84/4326 cause offset?

texts:
I am using the block of code below. I am reading in a geotif with a UTM zone projection and converting it to 4326. I am pulling out certain multiband values and creating a new raster, when I plot the new raster in ArcMap or QGis, the two don't line up. Can anyone tell me if I am doing this transform correctly?

Here is some additional information as requested
the offset is about 85 meters according to QGis
Geotools version is 11.1

ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY
            .createValue();
    policy.setValue(OverviewPolicy.IGNORE);

    // this will basically read 4 tiles worth of data at once from the disk...
    ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();

    // Setting read type: use JAI ImageRead (true) or ImageReaders read methods (false)
    ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
    useJaiRead.setValue(true);

    GridCoverage2D image = new GeoTiffReader(f).read(new GeneralParameterValue[]{policy, gridsize, useJaiRead});
    /**
     * reproject to WGS84
     */
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
    image = (GridCoverage2D) Operations.DEFAULT.resample(image, targetCRS);

    GridGeometry2D geometry = image.getGridGeometry();

    BufferedImage img = ImageIO.read(f);
    WritableRaster raster = img.getRaster();


texts: Watershed a catalog RasterStack in place of RasterLayer in lidR?

texts:
I want to use watershed with lidR package on a catalog. I paste this command:

rep <- "/home/lasfiles"
new_ctg <- catalog(rep)
lidR:::catalog_laxindex(new_ctg)
opt_output_files(new_ctg) <- paste0(rep, "/outputs/mnh/{ID}_cano")
opt_chunk_size(new_ctg) <- 500
opt_chunk_buffer(new_ctg) <- 30
opt_cores(new_ctg) <- 7L
algo <- pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
chm <- grid_canopy(new_ctg, 0.5, algo)

The result CHM is ok,

summary(chm)

        grid_canopy
 Min.       1305.953
 1st Qu.    1406.373
 Median     1468.121
 3rd Qu.    1539.043
 Max.       1716.370
 NA's        180.000
 Warning message:
 In .local(object, ...) :
     summary is an estimate based on a sample of 1e+05 cells (5% of
     all cells)

but when I pass this command:

algo <- watershed(chm, th = 4)

  Error in assert_is_all_of(chm, "RasterLayer") :
    chm is not a RasterLayer it is RasterStack.

Why this error?


texts: How to get layers from current viewport in Openlayers?

texts:
In my web map application, I have developed the printing module by using Mapfish with Geoserver. While printing, I need to send the layers to geoserver only which are falls inside the current view port.

Is it possible?


texts: Creating longitudinal profile for pipeline using QGIS?

texts:
I have a sewer system in QGIS and i wish to know if is possible to create a longitudinal profile between two manholes, showing the terrain and the pipeline.
In each point (manholes) I have in differents fields: ID, depth, Z terrain and each line has the diameter, but doesn`t matter if the pipeline is not represented with the real diamension and only as a line, the same for terrain, if is only a straight line between manholes is ok, cause the distances are too short and are in a city

Is it possible to create a layout with this graphic?


texts: LAStools errors in QGIS: 'cannot open lasreaderlas with file name .las', 'could not open lasreader'?

texts:
I am experiencing some problems with LAStools. I have two QGIS separate installations: 2.14.17 is a standalone installation and 2.18.11 is a Osgeo4w installation.

I have followed the procedure to use LAStools described here (in Spanish) and here. Each time I want to view a .las file (with lasview), the command seems to start but it closes (it happens with other tools as well).

In the log console I have the following message:

  017-07-27T08:37:10 0   LAStools command line
              C:\lidtools\bin\lasview -i "C:/Users/Mapoteca2/Desktop/demo/demo.las" -points 5000000
              LAStools console output
              ERROR: cannot open file 'demo.las'
              ERROR: cannot open lasreaderlas with file name 'demo.las'
              ERROR: could not open lasreader

This is odd since I am able to use the .las file from LAStools GUI. I have declared all LAStools executables as exceptions in my antivirus but the problem persists in both installations. I am running both installations of QGIS in Win10 64 bits with 32gb of ram.


texts: exportToPdf doesn't export the layout in QGIS 3

texts:
I have a trouble exporting a print layout in QGIS 3.10 to PDF when doing it programatically. I've created a project with print layout manually and now just want to modify the content of the map with python and then save it as PDF. This is the code I tried for export, everything works just fine except the last line, which does nothing and ends with code 4.
manager = QgsProject.instance().layoutManager()
layout = manager.layoutByName("my_layout")
exporter = QgsLayoutExporter(layout)
filename = "C://layout.pdf"
#so far so good, but this doesnt create the file:
exporter.exportToPdf(filename, QgsLayoutExporter.PdfExportSettings())

Any idea what could be wrong? When I run the export of the same layout manually, it works fine.


texts: Exporting results from Select By Location comparisons of each line in one shapefile with polygons in another?

texts:
I use ArcMap 10 and I am trying to export multiple selections (a.k.a 'Select by Location') from one file against another.

I start by selecting the first line of shapefile A and I want to select all the polygons of shapefile B that are inside of shapefile A (B polygons contained in A), to then export them as shapefile_1.

And this process to be repeated for x records in A to produce more shapefiles (shapefile_1, shapefile_2, ... shapefile_x).

How can this be done?


texts: How to create a point layer from a one-to-many join?

texts:
I have two layers:

One with geographical coordinates containing every address in the county and a unique identifier.
The second contains a list of incidents with the unique identifier from the address layer.

I could join the incidents layer to the address layer and then query only those addresses that have an incident. However there could be multiple incidents at on address, so I need to join the address layer to the incidents table. This gives me a record for each incident and the X Y coordinates from the address layer.

How do I specify that this layer now contains points without writing an intermediate file and re importing?

The two original layers will change with time and I don't want to save intermediate files.


texts: How to automate making thousands of raster layers from a CSV file in R?

texts:
I have a csv file of discrete point locations containing lat and long for different genera, with multiple entries for most genera, and I want to create a raster stack with a raster layer for each genus, with the rasters indicating the presence/absence of each genus. Because there are over 1700 genera in my file, I need to create some kind of loop to do this, unless there's an easier way, but for the life of me, I can't figure out how to do this. This is the code I wrote- I am aware that this is the wrong way to write loops, but I don't know what else to do. Here is my code:

cell_size <- 2.5
lon_min <- -179.0816; lon_max <- 178.5042; lat_min <- 24.5542; lat_max <- 81.817
ncols <- ((lon_max - lon_min)/cell_size)+1; nrows <- ((lat_max - lat_min)/cell_size)+1

for (i in cropaphidlocSPdata$Genus) {
  while !(i %in% aphidlist) {
    [i]data <- subset(cropaphidlocSPdata, Genus == [i]
    [i]count <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max,
                       res=2.5, crs="+proj=longlat +datum=WGS84")
    [i]coords <- cbind([i]data$Lon, [i]data$Lat)
    [i]count <- rasterize([i]coords, [i]count, fun = "count")
    aphidcounts <- addLayer(aphidcounts, [i]count)
    aphidlist <- c(aphidlist, [i])


texts: How to run a temporal query on a CSW using GeoNetwork?

texts:
On my GeoNetwork 2.10.3 a CSW GetRecords request (CQL) for TempExtent_begin after a certain date returns all of my metadata records. It seems like the filter expression is simply ignored.

http://myserver:8080/geonetwork/srv/eng/csw?SERVICE=CSW&version=2.0.2
&REQUEST=GetRecords&resultType=results
&constraintLanguage=CQL_TEXT&constraint_language_version=1.1.0
&constraint=TempExtent_begin%20after%202014-10-12T00:00:00Z
&elementSetName=full&outputSchema=http://www.isotc211.org/2005/gmd

Is this a problem of GeoNetwork or is there something wrong with my request?


texts: What's a good way store borehole data in SQLite / SpatiaLite?

texts:
I have a bunch of bores / wells with downhole data. The 'header' info for each borehole can easily go into a single table, but I'm wondering how best to store the downhole log data. Typically they are stored in spreadsheets as semi-continuous depth-value type data: e.g.:
depth,temperature
...values...

or there may be more than one value for each depth:
depth,temperature,gamma,neutron,...
...values...

I want to be able to easily pull them out to plot (e.g. XY scatterplots) in other software. Can I store entire text files in some sort of field?


texts:
Using GeoServer, I uploaded a netCDF file, got the temperature variable and mounted a color scale through the SLD.

This shows me a layer with the style I created.

That done, I'm presenting this layer in OpenLayers 3, but when applying a 'zoom in' everything is pixelated, without quality.

How can I fix the problem?

I noticed that in the GeoServer there is the bicubic interpolation option, but in doing so I can not pull the image on the map.

Below the excerpt of the relevant code:

//LAYER
var padrao = new ol.layer.Tile({
    source: new ol.source.OSM(),
    visible: true,
    name: 'padrao'
});

var temperatura = new ol.layer.Tile({
  extent: ext,
    visible: true,
  source: new ol.source.TileWMS({
    opacity: 0.5,
                url: "http://docker1.es.met.com.br/geoserver/wms",
                params: {'LAYERS' : 'somar:Temperatura',
                                 'TILED': true,
                             },
                serverType: 'geoserver',
                crossOrigin: 'canvas',
                projection: 'EPSG:4326',
            })
        });

        var view = new ol.View({
            center: [-6217890.205764902, -1710870.6048274133],
            zoom: 4,
            maxZoom: 7,
            minZoom: 4,
            extent: ext
        });

//Instancia o mapa com as configurações
var map = new ol.Map({
    target: 'mapa',
    controls: ol.control.defaults().extend([
        new ol.control.ScaleLine(),
        new ol.control.ZoomSlider()
    ]),
    renderer: 'canvas',
    layers: [padrao,temperatura],
    view: view,
});

I'm new with GeoServer and OpenLayers 3. Any direction of how to solve the problem will be of great help.


texts: Pyqgis equivalent of OID@ token (in arcpy cursors)

texts:
In ArcGIS I often find it very useful to use the OID@ (object id) token of for example da.SearchCursor when iterating over rows in vector data:

  Additional information can be accessed using tokens (such as OID@) in
  place of field names

For example if I for some reason would like to select last occurence of duplicate records I would in ArcGIS do (just as an example, made up data etc...):

d = {key:val for key,val in arcpy.da.SearchCursor("OIDtest",['sometext','OID@'])} #Dictionary cant have duplicate keys, only last record is kept.
d
{u'A': 1, u'B': 3}
d.values()
[1, 3]
arcpy.SelectLayerByAttribute_management(in_layer_or_view="OIDtest", where_clause="{0} IN{1}".format(arcpy.Describe("OIDtest").OIDFieldName, tuple(d.values())))

But in QGIS the OID column is hidden(?) and cant be accessed by name. Are there an OID@ equivalent that can be used to select row 2 and 4 when there is nothing else unique to identify the rows by?


texts:
I am using QGIS version 2.14.1 (Essen). I'm trying to do some spatial analysis, including an intersect between the same shapefile.

Basically I have a shapefile with the geometries of all the events planned in an area at a specific time. I'm trying to find the intersect between each geometry but when I do an intersect between the same shapefile, it returns all the geometries. I guess if I save every single geometry as a different shapefile then it is possible but as I have over 1000 data values I don't really want to do this.

Anyway I can do this in QGIS?


texts: How to insert layer to PostgreSQL from QGIS plugin

texts:
I need to insert the data of a layer of the qgis project into an existing table of postgis I have been using this code to test, but it gives an error since the method used creates the table in postgresql does not insert the elements if the layer exists.
def test_insertar_en_Postgis(self):
    app = QgsApplication.instance()
    self.conf = QgsSettings(
        app.qgisSettingsDirPath() + "QGIS/gestion_proyectos.ini", QSettings.IniFormat
    )

    self.conf.beginGroup("/conexion")
    self.servidor = self.conf.value("servidor")
    self.puerto = self.conf.value("puerto")
    self.conf.endGroup()
    self.conf.beginGroup("/cfg_nomenclatura")
    self.nomenclatura = self.conf.value("nomenclatura")
    self.directorio = self.conf.value("directoriogpk")

    self.usuario = self.cfgdlg.v_usuario.text()
    self.contrasena = self.cfgdlg.v_contrasena.text()

    layers = self.iface.mapCanvas().layers()
    indice =0;
    for layer in layers:
        indice= indice + 1
        con_string = "dbname='bdsapromat' host='" + self.servidor + "' port='" + self.puerto + "' user='" + self.usuario + "' password='" + self.contrasena + "' key=geo_id type= GEOMETRY  table= \"25_areas_verdes_terrenos\".test"+str(indice)+" (geom)"

        err = QgsVectorLayerExporter.exportLayer(layer, con_string, 'postgres',QgsCoordinateReferenceSystem('4326'),True)

        if err != 0:
            QgsMessageLog.logMessage('Hay error al insertar en postgis ', "Proyecto", level=Qgis.Info)


texts: File GDB Feature Classes Converting Themselves to File GDB Tables

texts:
I have GDBs in folders in a certain location in a shared dropbox. When I move the folders to another location, some of the feature classes have been converted to GDB tables.

I have never seen this before, and now all of a sudden it is happening to me all the time.

Has anyone else seen this and know what is causing it?


texts: How selecting uniques linestrings from two datasources with similar features?

texts:
I know is easy to select uniques points, or polygons features from two datasources (using centroids and overlap, or buffer) , but it not so easy do the same with linestrings, many times it is even a nightmare.

So,  how can you select uniques linestring from two datasources when there is very similar linestring features, some solaped, very near, o similar but not equals like this, or even more difficult, and if, in addition,  you also want to select the unique segments of those linestrings ?


texts:
I'm trying to create and if/else loop for three sublayers in a group layer. If it's the group layer (0) then the sublayers have certain labels. If it's group layer (4) then the sublayers have other labels, etc.

Here's what I've tried. I've also tried arrays [0] instead of '0'.

                   if(identifyParams.layerIds === '0'){

                          if(layerName === '10% Damage Risk'){
                               layerLabel = '24 Hour Forecast';
                               label = "10% Damage Risk";
                               }

                          else if(layerName === '90%  Damage Risk'){
                               layerLabel = '24 Hour Forecast';
                               label = "90% Damage Risk";
                               }

                          else if(layerName === 'Overall Damage Risk'){
                               layerLabel = '24 Hour Forecast';
                               label = "Possible Damage";
                               }

                          }

                    if(identifyParams.layerIds === '4'){

                          if(layerName === '10% Damage Risk'){
                               layerLabel = '48 Hour Forecast';
                               label = "10% Damage Risk";
                               }

                          else if(layerName === '90% Damage Risk'){
                               layerLabel = '48 Hour Forecast';
                               label = "90%  Damage Risk";
                               }

                          else if(layerName === 'Overall Damage Risk'){
                               layerLabel = '48 Hour Forecast';
                               label = "Possible Damage";
                               }

                          }

I've also tried, but it didn't work.

   if(layerName === '24 Hour Forecast')


texts:
I wrote a QGIS plugin including the "gdalogr:cliprastermymasklayer" algorithm but get the following error message:

That's the code where I used the algorithm:

    def clipRasters(self):
    inputDir = self.getInputPath()
    outputDir = self.getOutputPath()
    for lyr in glob.glob(os.path.join(inputDir, "*xyz")):
        processing.run("gdalogr:cliprastermymasklayer", lyr, self.getVectorLayer(), 'none', True, False, 0,
                       outputDir + lyr)

Though, the error message is not specifically refering to my usage of the algorithm, I guess that's where the problem is.

Now, when I look at the following QGIS Dokumentation https://docs.qgis.org/3.4/en/docs/user_manual/processing_algs/gdal/rasterextraction.html#clip-raster-by-mask-layer, it seems that I could even give 9 arguments. Thus, why does it say in the error message that I can use only 6? Or does the message not refer to the number of arguments for the gdal algorithm?


texts:
I´m new to Qgis and GRASS.
I need to change the default parameter skip of the module r.to.vect.lines, from its default value of 10 to 1.

How can I do this?

r.to.vect.lines [-v] input=name output=name [skip=integer] [--help] [--verbose] [--quiet] [--ui]


texts:
Does anyone know of a source for an Arctic EPSG:4326 map projection?

Obviously I'm a GIS newbie, so gentle pointers in the right direction appreciated!

It's a long story.  I have some arctic ice data that I'd like to visualize with Tableau 8.0 which accepts external wms maps, so long as they're EPSG:4326 (Plate Carree) projections.

Obviously such a projection tangent to the equator won't work for high latitudes.  But there's nothing that says you couldn't do a mercator-like projection using the prime meridan as the 'equator' so that close to the pole x-y placements would be reasonably accurate.

If such a thing existed I could upload it into Tableau 8.0 and plot my points.  Or, if there was a recipe for reprojecting an existing map into this format using grass or gdwalp or some such I could do it myself.

Anybody ever hear of such a thing?  Or have any advice?

Thanks.


texts: QVariant is not defined when using addAttributes

texts:
When I try to add a field programmatically with the following code:

import os.path
from PyQt4 import QtGui, QtCore
legend = qgis.utils.iface.legendInterface()
allLayers = legend.layers()
for i in allLayers:
  lyrName = i.name()
  GDBname = str(lyrName)
  print GDBname.split('_', 1)[0]
  lyrProv = i.dataProvider()
  i.startEditing()
  lyrProv.addAttributes( [ QgsField("testName", QVariant.double) ] )

I receive an error:

lyrProv.addAttributes( [ QgsField("testName", QVariant.double) ] )
NameError: name 'QVariant' is not defined

Can someone help?


texts: finding the correct transformation (from EPSG:31300 to EPSG:4326) or how to store two srid in one table

texts:
I have a very specific problem. I have data crossing two coordinate-systems, more specifically: the Netherlands and Belgium.

Previously the data was stored in one table in Oracle, using the RD coordinate system. When showing the data from Belgium, it was converted using the Oracle command SDO_CS.TRANSFORM back to the Belgian coordinate system (Lambert 72) and all was well.

Now I am moving to geoserver, and rendering using openlayers, on top of openstreetmap or google maps background. I notice the Belgian data (which is stored in RD) is shifted. If I create a table with the original Belgian data (in Lambert 72 without transformation), it is rendered correctly.

So my first guess: the Oracle transformation from LB72->RD is wrong.

But if geoserver does the transformation it is correct (since it renders the original Lambert 72 data correctly on the google maps background).

So I imported the data from one table to Postgis, and do the transformation there (using ST_transform(geom, 4326) and I see it is even more wrong :)

So at this moment I have a problem:

I have data in two coordinate systems, which belong in one table since it is border-spanning (it is a pipeline)
the data is rendered correctly, by geoserver, if using the original coordinate system
it is not if the belgian data is "transformed" using database tools (SDO_CS.TRANSFORM or ST_TRANSFORM)

My question:

how can i correctly transform belgian data to a coordinate system that can be rendered correctly on top of google maps or openlayers (I was considering WGS84 or ETRS89, as a shared coordinate system for both my RD and Lambert72 data)
or alternatively: how can I store two coordinate-systems in one table? (it is not possible in Oracle, is it possible in Postgis?)
or alternatively: how does one generally handle data that spans several coordinate systems?

Any help is greatly appreciated :)


texts: Recommended instance type when using pgrouting in AWS

texts:
Some advice required please.

When choosing instances in AWS (amazon web services) I can have a 'compute optimised' or 'memory optimised' instance. Can anyone recommend the best to use with pgrouting?

I am using pgrouting 2.4.2 with PostgreSQL 9.6.6


texts: Getting origin key to automatically populate foreign key while using relationship class in ArcGIS for Desktop?

texts:
I am using ArcGIS for Desktop, and I am attempting to create feature linked annotations for my cable system and I have set up a relationship between the two features so that any changes made to the cables will automatically be updated in the annotation.

The cable feature's "GUID" field is set as the origin and I have created a "CableGUID" field in the annotation feature acting as the foreign key.

I can select a cable and choose to "Annotate Selected Features" which will generated the annotation but the origin key is not being migrated to the foreign key column.

Am I missing something with the set up?


texts: IFeature.Store seems to bypass edit sessions and cannot be rolled back. Why?

texts:
I'm working with ArcObjects edit sessions for the first time and am having some trouble with it. The C# code I am about to discuss is part of an ArcMap 9.3 extension DLL.

My problem is that calls to IFeature.Store inside edit operations seem to simply bypass the edit operation (and edit session) in which they execute.

This is my code:

IFeatureWorkspace workspace = …;
IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)workspace;
if (!workspaceEdit.IsBeingEdited()) return;

// at this point we know that we are inside an edit session.

IFeatureClass featureClass = workspace.OpenFeatureClass(…);

workspaceEdit.StartEditOperation();
try
{
    IFeature feature;
    IFeatureCursor features = featureClass.Search(null, false);
    // (ESRI recommends search cursors inside ArcMap, even when updating inside edit sessions.)
    while ((feature = features.NextFeature()) != null)
    {
        feature.set_Value(…); // only as an example; assume that the feature is modified somehow.
        feature.Store();
    }
    workspaceEdit.StopEditOperation();
}
catch
{
    workspaceEdit.AbortEditOperation();
    throw;
}

As you can see, the edit operation in this code modifies all features in some feature class.

Let's say that after a few successful feature updates, one update throws an exception. The exception handler in the above code would then call workspaceEdit.AbordEditOperation();. As I understand it, this should rollback all modifications done inside the edit operation (including the previous, successful feature updates). I am however observing that the effects of feature.Store(); never get rolled back.

I have also tried this with IEditor.StartOperation/.StopOperation/.AbortOperation instead of the IWorkspaceEdit methods, but that doesn't change the observed behaviour.

Can anyone explain to me what I am doing wrong? How can I rollback the effects of IFeature.Store in an edit session?


texts: Creating classification matrix for rasters with abrupt class intervals in R

texts:
I am trying to re-classify to raster images which do no have continuous range intervals which normally have been used in examples for creating classification matrix. My original data has these values and classes it has 17 classes in total and i want to aggregate and reduce them for better visualization of data:

Needleaf Forest 1
Broadleaf Forest 2
Deciduous needleaf Forest 3
Deciduous broadleaf forest 4
mixed Forests 5
closed shrublands 6
open shrubland 7
woody savannas 8
savannas 9
grasslands 10
wetlands 11
croplands 12
built up 13
crop/natural veg 14
permanent snow 15
barren 16
water bodies 17
unclassified 255

Now I want to reclassify this according to following table:

1-5-> 1
6-7->2
8-9->3
10->4
11->5
12,14->6
13->7
15,17->8
16->9
255->NA

I am using the following code to built up the matrix and reclassify images in batch mode but the matrix is not able to classify when it's not given in a sequential order. How can i produce the above mentioned classification?

Here is the code i am working on:

library("raster")
grids <- list.files("/media/MyData/Data/MODIS/PROCESSED/MCD12Q1.006_20190728202245", pattern = "*.tif$")

reclass_df <- c(1, 5, 1,
                6, 7, 2,
                8, 9, 3,
                4,255, 10)
rclmat <- matrix(reclass_df, ncol=3, byrow=TRUE)
r <- raster("/media/MyData/Data/MODIS/PROCESSED/MCD12Q1.006_20190728202245/MCD12Q1.A2003001.LC_Type1.tif")
r
rc <- reclassify(r, rclmat)
rc
#function to reclassify rasters and write a new reclassified tif file for each
batch_reclass <- function(grids){
  for (i in 1:length(grids)) {
    #read in raster
    r <- raster(paste0("/media/MyData/Data/MODIS/PROCESSED/MCD12Q1.006_20190728202245/", grids[i]))
    #perform the reclassifcation
    rc <- reclassify(r, rclmat)
    rc
    #write each reclass to a new file
    plot(rc)
    writeRaster(rc, filename = paste0("/media/MyData/Data/MODIS/PROCESSED/reclass/", "rc_", grids[i]), format=".tif")
  }
}
#run the function
batch_reclass(grids)

The current matrix in code works fine but something like:

c(1,5,1,
                6,7,2,
                8,9,3,
                10,4,
                11,5,
                12,14,6,
                13,7,
                15,17,8,
                16,9,
                9,255,NA)

It throws up error

Error in .local(x, rcl, ...) : all(rcl[, 2] >= rcl[, 1]) is not TRUE

which is because of incorrect matrix. How can i create the right matrix?


texts: Display Leaflet Markers with UTM coordinates?

texts:
is it possible to display markers with UTM coordinates (zone 32) with leaflet and OSM? How can I do this? If I use UTM Coordinates like latlng, the position is totally wrong:

        L.marker([5698015.694, 312400.093]).addTo(map)
        .bindPopup("Hello world!
I am a popup.").openPopup();

Did I have to translate them first, or is it possible to display UTM coordinates directly with leaflet?


texts: Cannot edit nodes of layer imported from KML

texts:
I'm new to QGIS. This group has helped me understand how to import KML files as vector layers, but once the data is imported, I cannot figure out how to edit the nodes of the imported polygons.  I can see them as red Xs when I toggle edit mode, but the node tool is greyed out/disabled.

How can I edit the nodes of these imported polygons?


texts: How do I run the nvis tool in QGIS GRASS?

texts:
I am trying to follow this workflow: [How to generate a DEM from a contour Shapefile?] and it's running fine, and when I hit view output or view results I end up with a flat raster with shading associated with the heights of my contours on screen.

However, I want to see the 3D terrain with hillshading etc, and potentially export a 3d mesh, and from research, think that the nvis module in GRASS can do this for me, but when I try to run this I get the error "module nvis not found"? And QGIS crashes and dies.

When I go to Manage and Install Plugins, Grass is there, and it has let me perform all the other modules I've needed? So not sure what to do with this??

Any ideas or help much appreciated!

Edit: Link to sample data here (.shp+.gml): https://www.dropbox.com/sh/mwxlr3ij78d8gkr/lKwgu3Pafu


texts:
I am a novice to Python and QGIS coding. I am facing an issue while showing the MySQL spatial data on QGIS canvas. Database connection details are below:

database name:- mysql
host :- localhost
port :- 3306
user :- root
password :- abcx
tablename:- Tablegis

I found a code on to connect the same:

uri="MySQL:mysql,host=localhost,port=3306,user=root,password=abcx|layername=Tablegis"
vlayer=QGSVectorLayer(uri, "Open_Table","ogr")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

This code doesn't show any warning or error but nothing shows on my canvas. When I am printing the vlayer shows the object reference. Tried hard nothing worked for me.


texts: Arcpy error: ERROR 000732: Clip Features: Dataset xxx does not exist or is not supported

texts:
I've been working on a project. A small part of it is to clip point features by ellipse polygons.

I need to loop through a set of IDs, when I call "clip" function via arcpy, a very strange error persists:

"ERROR 000732: Clip Features: Dataset Ellipse/ellipse_6406685702.0 does not exist or is not supported
Failed to execute (Clip)"

the clip function worked fine if only a single ID was used, the error only occurs when it loops through the IDs.

I've checked Esri's help, and tried to appy some fixes, but none of them  worked. I'm pretty sure the files such as ellipse_6406685702.0 indeed exist as I can use them one by one, but cannot use them in the batch way. This is the very last step of my project so I really need to get out of it.

Here is the code:

import arcpy
arcpy.env.workspace= "C:/Users/Ran/Desktop/TucsonTravelAnalysis/data/TucsonTravel.gdb"
fc = "All_Locations/Tucson_ppl"
pplWH= []
cursor = arcpy.da.SearchCursor(fc,["ID_individ","TripType"])
for row in cursor:
        if row[1]== 2:
            pplWH.append(row[0])
pplWH= list(set(pplWH))
for ID in pplWH:
 infc = "All_Locations/All_ND_points"
 clipfc= "Ellipse/ellipse_{0}".format(ID)
 outclipfc= "All_Locations/HW_NDpoints_{0}".format(ID)
 arcpy.Clip_analysis(infc, clipfc, outclipfc)


texts: How to query all points within the current map extent

texts:
First post - apologies if it's a simple question!

I am basically looking to have a JS API service which will populate a datagrid with all of the points visible in the current map extent. As a user zooms in or out, the number of records returned will change to reflect what is in the current extent.

I have been looking at the findPointsInExtent examples in the esri JavaScript examples, but am not sure how to edit the these examples to go from a user defined rectangle to the map extent as described above.

Many thanks,


texts:
I try to display the France map on HTML with RStudio; I followed the tutorial to the World card in Creating an interactive map with Leaflet; Example works.

I wanted to do the same just by displaying the map of France:
code:

library(rgdal)
library(foreign)
library(maptools)

FRA=readShapePoly("/Users/insa/Documents/exempleR/FRA_adm/FRA_adm0.shp")

#plot(FRA_adm0)
# Nettoyage des données

library(leaflet)
## Initialisation
m <- leaflet(padding = 0)
## Ajout des pays

## Ajout des pays
m <- addPolygons(map = m, data = FRA, opacity = 100,
                 color = "#FAFCFA",
                 weight = 0.25,popup = NULL,
                 options = list(clickable = FALSE),
                 fill = T, fillColor = "#B3C4B3",
                 fillOpacity = 100)

## Dimensions de la carte
m$width <- 874
m$height <- 700

# Export de la carte en html
library(htmlwidgets)
saveWidget(m, 'mapfrance.html', selfcontained = TRUE)

error:

Stack space overflow: current size 33632 bytes.
Use `+RTS -Ksize -RTS' to increase it.
Erreur : pandoc document conversion failed with error 2

I do not understand this error, is that it's the input file is shapefile .shpaand not .rdata?


texts: Clipping a polygon shp by one large polygon

texts:
I've seen similar answers to this question on so many forums but nothing seems to work.

I have one large polygon (a country boundary) with another polygon layer that contains thousands of smaller polygons (separate sites, either independent polygons or sharing the same boundary) that fall inside, outside and across the large polygon.
I want only the small polygons that fall outside the larger polygon even if it breaks the smaller polygon that fall across the larger polygon boundary. I have used the geoprocesing->clilp tool but I just get an error message (i.e. the ! symbol) and no new layer is created.

I'm using Arc 10.3.1.
Is there any other tool I can use? (I've also tried intersection and select by location but clipping seems to be the one to use and it doesn't work!)


texts: Problem with QgsMapCanvas().zoomToSelection() when running from function

texts:
I have the following function which, when typed in to the python console line for line works fine, but when I use the function it gets as far as highlighting/selecting the accident (feature) but the zoomToSelection() statement doesn't appear to do anything.

def find_acc(layer, expression, scale = 2500):

    lay = core.QgsMapLayerRegistry.instance().mapLayersByName(layer)[0]
    canvas = QgsMapCanvas()

    request = core.QgsFeatureRequest().setFilterExpression(expression)
    iter = lay.getFeatures(request)

    selection = [f for f in iter]

    if len(selection) == 1:
        print "Found accident"
        accident = selection[0]
    else:
        raise IOError, "More than 1 entry found"

    lay.setSelectedFeatures([accident.id()])
    canvas.zoomtoSelected()
    #canvas.zoomScale(scale) disabled to see if would make a difference - it didn't

    # added in afterwards
    canvas.refresh()

once imported in to QGIS I use the function find_acc("accidents", "\"Police_ref\" = '0111XH30035'", 2500).

I tried again with canvas.refresh() after zoomToSelected() but that didn't help. Interestingly, despite the feature being highlighted on the map I can't use the keyboard shortcut CTRL + J or View > Zoom to Selection either.


texts:
I am currently having a problem with my image on OSM via OpenLayers. I can present my Image on the right spot (nearly exactly there where i want it).
The problem is, while zooming in or out, the image changes it's location, and how it should look like. On the zoomlevel I'm setting, the image looks like it has to.

Here is the code which i am using, to show my image on the right location:

var layers = [
        new ol.layer.Tile({
            source: new ol.source.MapQuest({layer: 'osm'})
        }),
        new ol.layer.Image({
            extent: [1447153.38031, 5780349.22026, 2003750.83428, 6446275.84102],
            source: new ol.source.ImageWMS({
                url: 'file:///C:/Users/Xhy/Desktop/linear.png',
                params: {'LAYERS': 'modis,global_mosaic'}
            })
        })
    ];

With the extent i bound my image on the coordinates I need. I had to transform my coordinates from WSG84 to Google Transverse Mercator, for that i used this site: http://cs2cs.mygeodata.eu/

The following code, is how my layers are shown on the map:

var map = new ol.Map({
        layers: layers,
        target: 'map',
        view: new ol.View({
          center: [1728012.45558,5982976.33295],
          zoom: 7
        })
      });

The center is also in Google Transverse Mercator, and i also set the zoomlevel on 7.

The output i get is nearly that what i want:

When I zoom in it changes:

When i zoom in once more, the change is drastic, the image is completely somewhere else:

So my question is, how can i bound my image on the right location without it changing while zooming in? Do i have to set the bound different for every zoomlevel?

It took me quite some time to get the image where I wanted, now I struggle with this problem, can anyone give me any tips on how to solve this?


texts: Create filter expression by selected elements?

texts: Is it possible to automatically create a filter expression for the Query Builder by selected elements?


texts: OpenLayers does not appear in the QGIS plugin menu

texts:
I have installed the OpenLayers QGIS plugin, which features external map layers in QGIS. Quoting the official page:

  OpenStreetMap, Google Maps, Bing Maps layers and more

In QGIS, it seems that the plugin is well installed:

However, no new items were added to the plugin menu:

Any idea how to start the plugin? I have tried restarting QGIS and reinstalling the plugin several times, with no success.


texts: Incremental "auto save" in ArcMap edit session

texts:
Is there any way to add automatic incremental saving to ArcMap 10.1 editing sessions?  I found this script for 9.x but it doesn't seem to work with the newer versions.  I'm looking for a script or add on. I'm more concerned with saving edits rather than saving the mxd itself.

Frequent crashes have ruined my morning too often.

I'm aware that this is similar to Does ArcGIS for Desktop have auto-backup feature?, but it is not exactly the same in my opinion.


texts: Is there more than one option to change QGIS' temporary folder?

texts: I have tried changing the TEMP folder of QGIS following this post,and although it wasnt required, I restarted the system after making the change, but still got no result (printscreen below). Is there an alternative way to change the location of the temporary folder for QGIS? Perhaps something about altering some installation files?


texts: .DWG and .DGN will not drop in to ArcMap 10.2.2

texts:
I am attempting to bring in Microstation data (.DGN) into a working ArcMap file. Every time I attempt to drop the data in the map I get the following error message... (Error could not open the specified file)

The original files were .DGN, I had a friend convert them to .DWG and have received the same error. I have also attempted to use the {CAD to Geodatabase} tool on both file extension types and, once again, get the same error message when I try to add the data to the map.

Any suggestions?


texts: How to implement VisibleBoundsUpdatedEventHandler

texts:
I have an ArcObjects add-in button that adds attribute value info as a text graphic into a map document layout.  Now, instead of using the add-in button click event, I would like the code to run on VisibleBoundsUpdatedEventHandler or when the active view has been updated.  Can this be implemented through add-ins?  The link below uses this for a panning tool.  I'm not sure how and where to begin to implement this.  Any suggestions or code samples would be most appreciated!

http://forums.esri.com/Thread.asp?c=93&f=993&t=145041

Thank you


texts: Find closest target in a network using GraphHopper

texts:
A person is somewhere in a city and wants to know the network distance to the closest supermarket.

Is there any way I could solve this problem with the GraphHopper program/library (not the API)?


texts: Changing text of ArcGIS Pro text element between repeated exportToPDF of layout?

texts:
I am using ArcGIS Pro 1.3.1 and it seems like ArcPy is not always respecting the current setting of the text property for a text element when I export its layout to PDF.

To reproduce this I have created the following test:

Start ArcGIS Pro and create a New project using the Blank template in C:\Temp called TestProject
INSERT a New layout using A5 from the ISO - Landscape category
INSERT a Text element somewhere in the middle of the page.
Leave the text as "Text" but give it a size of 72pt if you want to make it very easy to read.
Save the project and close ArcGIS Pro
Run the Python script below from the IDE of your choice

import arcpy

aprx = arcpy.mp.ArcGISProject("C:\Temp\TestProject\TestProject.aprx")
lyt = aprx.listLayouts("Layout")[0]

for elm in lyt.listElements("TEXT_ELEMENT"):
    elm.text = "TEST1"
lyt.exportToPDF(r"C:\Temp\TEST1.pdf")

for elm in lyt.listElements("TEXT_ELEMENT"):
    elm.text = "TEST2"
lyt.exportToPDF(r"C:\Temp\TEST2.pdf")

I would expect that this would create two PDF files with the text "TEST1" written in the one called TEST1.pdf and "TEST2" written in the one called TEST2.pdf.

Instead what I get is:

"TEST1" in TEST1.pdf (as expected)
"TEST1" in TEST2.pdf (I suspect that this is a bug)

It seems like the layout object is only allowing this text element to have its text changed once between repeated exports of a layout to PDF.

At the moment I am unable to run the same test using ArcGIS Pro 1.4 (where I am hoping it is fixed), and so I am wondering whether anyone can run it on their version of ArcGIS Pro to confirm whether they can reproduce it?

It seems to be the exportToPDF that "locks" the text value because if the second test below is run "TEST3" gets exported to TEST3.pdf OK, but TEST4.pdf gets "TEST3" instead of "TEST4".

import arcpy

aprx = arcpy.mp.ArcGISProject("C:\Temp\TestProject\TestProject.aprx")
lyt = aprx.listLayouts("Layout")[0]

for elm in lyt.listElements("TEXT_ELEMENT"):
    elm.text = "TEST1"
for elm in lyt.listElements("TEXT_ELEMENT"):
    elm.text = "TEST2"
for elm in lyt.listElements("TEXT_ELEMENT"):
    elm.text = "TEST3"
lyt.exportToPDF(r"C:\Temp\TEST3.pdf")
for elm in lyt.listElements("TEXT_ELEMENT"):
    elm.text = "TEST4"
lyt.exportToPDF(r"C:\Temp\TEST4.pdf")

As a sanity check I ran the ArcMap equivalent test below using ArcGIS Desktop 10.4.1, and that behaved as expected:

import arcpy

mxd = arcpy.mapping.MapDocument(r"C:\Temp\test.mxd")

for elm in arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT"):
    elm.text = "TEST1"
arcpy.mapping.ExportToPDF(mxd,r"C:\Temp\TEST1.pdf")

for elm in arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT"):
    elm.text = "TEST2"
arcpy.mapping.ExportToPDF(mxd,r"C:\Temp\TEST2.pdf")

Consequently, the ArcGIS Pro 1.3.1 behaviour mentioned at the beginning is looking more and more like a bug.


texts: Documenting custom Python functions in QGIS

texts:
How can I document my custom Python functions in QGIS?

I would like to offer users help for functions I wrote for them:

I've tried adding a docstring:

But whenever I click Load an error message is displayed:

The error message clearly shows an indentation error, which does not correspond to my actual code.


texts: Finding min and max symbology for each band at command line using GDAL?

texts:
I would like to automate the process of getting the correct symbology min/max scale for each of the RGB bands via command line.

My use case is I am converting multi-band GEOTIFFs to JPG using the following command

 gdal_translate -ot Byte -of JPG -b 1 -b 2 -b 3 -scale_1 0 439 -scale_2 0 460 -scale_3 0 454 -a_nodata none test.tiff output2.png -co PHOTOMETRIC=RGB

The scales change frequently, and without the correct min and max, the images appear black.

My solution (which does not scale) is I have been getting the scales by loading my GEOTIFFs into the QGIS3 UI software and then looking at what QGIS3 sets for the symbology mins and maxs for each color band. I would like to automate the process of figuring out the color band min/max scales via command line.


texts: Getting halfway point in Polyline in ArcGIS API for JavaScript?

texts:
I would like to reach the mid point of a polyline.

I just tried the centroid and realcentroid properties but I only got a soluction in ArcPy with the code below:

Midpoint = geometry.positionAlongLine(0.50,True).firstPoint
... print Midpoint.X
... print Midpoint.Y

How do I find the correct mid point (not centerline or gravity point of an irregular polyline)?


texts: How to disconnect all open GDB locks in ArcMap?

texts:
I am having trouble releasing GDB locks in ArcMap 9.3.

Does anyone know a way to disconnect from all open GDB connections in ArcMap 9.3.

I could close and re-open arcMap, but that seems to be a little time consuming.

Has anyone made a button to do this?


texts:
I create a tiff file from a bufferedImage in Java. I have some GCP (Ground Control Points) and I would like to make a GeoTiff. I know the GeoTIFF jai but this library is not updated since 2001...
Any ideas ?


texts: SQL Server error when creating Geography data from a Multipolygon Geometry

texts:
So, in SQL server, I have a table called Municipios, which contains a column with 2463 entries of geometry type, with their SRID set to 4326.
I'm running the following Query to try to convert them to geography type:
SELECT TOP 1 mun.nom_mun, geography::STGeomFromText(geom.STAsText(),4326) AS geog
FROM [dbo].[Municipios] AS mun

I also tried the following:
SELECT TOP 1 mun.nom_mun, geography::STGeomFromText(geom.MakeValid().STAsText(),4326) AS geog
FROM [dbo].[Municipios] AS mun

The problem with both is that the generated geography is definitely wrong, here are the original geometry and the resulting geography of the first entry in Municipios:
Geometry:

Geography:

I also, tried getting the string of the geometry and using that to build the geography, using:
SELECT mun.nom_mun, geom.MakeValid().STAsText() AS geog
FROM [dbo].[Municipios] AS mun

to obtain the Multipolygon string and:
SELECT TOP 1 mun.nom_mun, geography::STGeomFromText('MULTIPOLYGON (((...)))',4326) AS geog
FROM [dbo].[Municipios] AS mun

to generate the geography (I'd put the whole MULTIPOLYGON (((...))) here, but it's extremely long).
But, i still get the same wrong result.
I tried to cheat the system by using the Multipolygon string in Postgis (pgAdmin) and then move the generated geography back to SQL server. Inside pgAdmin, the resulting geography was correct:

But when i moved it back to SQL server, it was the same wrong geography.
Ideally, I'm looking for a way to do this conversion inside SQL server, if not, find a way to do it in postGis and move it back to SQL server without it becoming the same wrong projection, if none of those are options, well, any other way is good enough.


texts: Capturing geotagged text notes for field work?

texts:
I'm looking for a smartphone application that will allow me to take notes at a site in the field and associate each note with a geotag (like photos) of the current location which I can then later import in to GIS software.

What options do I have?


texts: Find a centerline of an intersection between two irregular shaped polygons

texts:
I have two irregular shaped polygons that overlap each other at least in at two points. I want to first find the intersection area between all of them, then for the polygon pair that intersects I want to find the centerline of each intersection, which should not be a straight line. I thought of using the Centerline module in python that uses Voronoi daigram, but it doesn't seem to work well. Is there another library that can quickly calculate the centerline of an intersection that starts at one intersection point and ends at another?

I have the following code:

import geopandas as gdp
from shapely.geometry import Polygon
from centerline.geometry import Centerline

crs = {'init': 'epsg:4326'}
attributes = {"id": 1, "name": "polygon", "valid": True}

# find the intersection of polygons
intersection_list = []
for x in range (0, len(list_of_polygons)):
    if len(list_of_polygons) > 1:
          # check if the last polygon in the list intersect with the current polygon
         if list_of_polygons[-1].intersects(list_of_polygons[x]):
                        polygon_last = gdp.GeoDataFrame(index=[0], crs=crs, geometry=[list_of_polygons[-1]])
                        polygon_current = gdp.GeoDataFrame(index=[0], crs=crs, geometry=[list_of_polygons[x]])
                        # get the intersection geometry
                        intersection = gdp.overlay(polygon_last, polygon_current, how = "intersection")
                        intersection_geom = Polygon(intersection.geometry)
                        # find the centreline of the intersction
                        centreline = Centerline(intersection_geom, **attributes)


texts:
I am using ArcGIS Pro for my map automation. I try to use arcpy.mp module to identify the maps in the current project. I want to export the map into a jpeg image like what we can do via Map Export under Share tab. But I didn't find any relevant methods under map object and as to my research, either layout or MapFrame could have method exportToJPEG. below is my code so far and no layout in my current project.

How can I use Python to convert a map object to MapFrame or layout object for exporting as an image or any other way to achieve my purpose?

aprx =
arcpy.mp.ArcGISProject(r"C:\Users\Documents\ArcGIS\Projects\myProject.aprx")

mps = aprx.listMaps()

lyr = mps[0].listLayers()

aprx.listLayouts()


texts: Creating a Webapp widget that pulls data from underlying polygons

texts:
I am currently wondering what would be involved in creating a custom widget that does the following.

Widget would be based off a point feature class. When you create a new point the information from the underlying polygon (specific fields chosen) will be appended to a field in the point feature dataset. e.g. If I have parcel fabric (polygon) displayed and i draw a point on a property the newly created point would inherit the address (attribute info) from parcel fabric (polygon).

An excellent example of what I am trying to accomplish can be viewed in the workforce for ArcGIS environment.

If you modify the assignments map that the application works off of and add a dataset (e.g. parcels or census, etc.) when you create a new assignment (point feature class) the pick point will append the popup info data from the polygon to the description field of the assignment dataset (point feature class).

Alternatively the geo-form app would be ideal, the only difference is when you select your location it would pull the polygon data into your point feature attributes as well as let you fill in your own fields as well.

I realize that what I am doing may not be possible. I figured this would be the best place to get the conversation started.


texts: Convert polygon co-ordinates to GeoJSON using JavaScript


texts: Can I query in ArcGIS Online for features that have attachments

texts:
I have a feature service in ArcGIS Online that has attachments enabled.  We upload photos of sidewalk trip hazards that have been painted by the city.

I have several queries setup with the Query widget to search for the date sidewalks were painted and other filters, but I want to be able to search for features that have photos attached.

Is there a way to query for and quickly locate features that have attachments without having to click through them until I see one with an attachment link?


texts: Adjust altitude of points in kml to ground level

texts:
I have recorded a path using GPS and imported the kml file into google earth.
Due to inaccuracies in recording the path (I assume) the points are sometimes above or below the ground.
Is there a way to adjust the altitude of the points to ground level?
In the path properties, I see the option to "Extend to ground", but this is not what I'm looking for.


texts:
I am trying to integrate an attribute query for an L.esri.featurelayer I'm currently using in one of my maps.

Currently I'm using this but it doesn't seem to work. The panel appears on the map but I can't seem to activate the selection of the various project types I'm trying to map out.

// Query to select project types
var ProjectTyp = document.getElementById('ProjectTyp');

ProjectTyp.addEventListener('change', function(){
stops.setWhere(ProjectTyp.type);
});


texts: Rapideye 3A Dark Object Subtraction in QGIS?

texts: Is there a method to do Dark Object Subtraction (DOS) with the Semi automatic classification plug-in of QGIS on Rapideye 3A imagery?


texts: Using rotated projection in NetCDF file?

texts:
Could anyone show me a working example of a NetCDF file header with a rotated Cartesian projection where the datum is WGS84?

From Spatial reference for netCDF data  I see that coordinates can be called  projection_x_coordinate and projection_x_coordinate but it doesn't specify how to fill in the projection type.

Looking at CF Standard Name Table I see variables such as "angle_of_rotation_from_east_to_x" but I don't know the full set necessary to define a rotated projection.


texts:
The following python script gives me a choropleth of the maternal mortality rate by country.

The resulting map, though, assigns the same color to (e.g.) Congo and Canada, which is not what the data say. I'm trying to learn how to do this, and I'm stymied by this strange result.

import requests
import zipfile
import io
import os
import pandas as pd
import geopandas as gpd

mmr = requests.get('http://api.worldbank.org/v2/en/indicator/SH.STA.MMRT?downloadformat=csv')
mmr_zip = zipfile.ZipFile(io.BytesIO(mmr.content))
names = zipfile.ZipFile.namelist(mmr_zip)
mmr_csv = zipfile.ZipFile.extract(mmr_zip, names[1])
mmr = pd.read_csv(mmr_csv, header =  2)
mmr = mmr[['Country Name', 'Country Code', '2015']]
mmr = mmr.rename(index = str, columns = {'2015': 'mmr'})

# get world shapefile with country borders
world = requests.get('http://thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip')
world_zip = zipfile.ZipFile(io.BytesIO(world.content))
zipfile.ZipFile.namelist(world_zip)
world_zip.extractall()
world_map = gpd.read_file('TM_WORLD_BORDERS-0.3.shp')

merged_map_mmr = pd.merge(world_map, mmr, left_on = 'ISO3', right_on = 'Country Code')
merged_map_mmr.plot(column = merged_map_mmr['mmr'], cmap = 'RdYlBu', scheme = 'QUANTILES', figsize = (10, 8))


texts: Using ST_MakeEnvelope over the antimeridian

texts:
I'm using ST_MakeEnvelope in PostGIS to try and make geometries that cross the international dateline (antimeridian), but they are making a shape that wraps the other way around the globe.

The geometry I want is one which covers the eastern half of Australia and the dateline.

It looks something like:

ST_MakeEnvelope(133, -49, -170, 0, 4283)

But:

ST_Intersects(ST_SetSRID(ST_MakePoint(150, -25), 4283), ST_MakeEnvelope(133, -49, -170, 0, 4283));

returns false. The shape goes the wrong way around the world!

If I just take:

ST_Intersects(ST_SetSRID(ST_MakePoint(150, -25), 4283), ST_MakeEnvelope(133, -49, 180,v0, 4283));

I get true, as expected.

How can I set a direction to my envelope? Or should I create a polygon some other way?

(If you're wondering: I'm trying to select geometries in a table for a map view, and when the view includes the dateline, nothing is shown on the map)

[ETA]

If it helps, the box actually arrives in Spherical Mercator, looking something like:

14855319.859171, -6294627.0589708, 21149946.918142, 0

But calling ST_Transform:

st_transform(st_makeenvelope(14855319.859171, -6294627.0589708, 21149946.918142, 0, 900913), 4283)

Results in geometries with the negative-valued longitudes:

SRID=4283;POLYGON((133.447608799776 -49.1163523261112,133.447608799776 0,-170.006794251165 0,-170.006794251165 -49.1163523261112,133.447608799776 -49.1163523261112))


texts: projectRaster() result does not have the same nrow and ncol, but has same resolution

texts:
I have a land use land cover (lulc) raster in WGS84 coordinate system shown below

class       : RasterLayer
dimensions  : 59575, 57541, 3428005075  (nrow, ncol, ncell)
resolution  : 0.000509, 0.000509  (x, y)
extent      : 68.12366, 97.41203, 6.754034, 37.07771  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source :<address path>
names       : lulc raster
values      : 1, 19  (min, max)

I have another raster in MODIS sinusoidal projection below which I will refer to as the foundation raster

class       : RasterLayer
dimensions  : 7276, 8330, 60609080  (nrow, ncol, ncell)
resolution  : 463.3127, 463.3127  (x, y)
extent      : 6508617, 10368012, 751956.5, 4123020  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs

data source :
   names       : foundation raster
   values      : 2, 3  (min, max)

I have to do analysis with both rasters and cannot change the projection of the foundation raster. So I reprojected the lulc raster as below

lulc_projected<-projectRaster(lulc raster, res=res(foundation raster),crs=crs(foundation raster), method='ngb') #categorical classes; hence ngb.

The resulting raster lulc_projected is completed without errors and aligns perfectly with the foundation raster when viewed in QGIS and Arc as below.

 class       : RasterLayer
 dimensions  : 7258, 10179, 73879182  (nrow, ncol, ncell)
 resolution  : 463.3127, 463.3127  (x, y)
 extent      : 6055308, 10771368, 744727.8, 4107452  (xmin, xmax, ymin, ymax)
 coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : lulc projected
values      : 1, 19  (min, max)

However, the nrows and ncols differ not allowing me to continue with analysis such as mask() in R. More confusing is that the resolution matches.
I know I can change the extent by using extent(lulc_projected)<-extent(foundation raster).

As a solution I resampled the lulc_projected raster as below

lulc_projected_resampled<-resample(lulc_projected,foundation_raster, method='ngb')

After this the rasters match perfectly i.e. extent, resolution, coordinate system and nrows and ncols. But when I view both the lulc_projected_resampled and foundation raster, there is a shift by about 40 km. Refer to this post for a visual example of this shift and comments (with no solution).

So, firstly, why is projectRaster() not giving a raster with the same number of rows and columns? The product at this stage i.e. lulc_projected aligns perfectly with the foundation raster, however I cannot proceed with more analysis in R, such as stack() and mask() because the nrows and ncols do not match.

Secondly, is this not the workflow? If not, what is the workflow to project a raster to match another?


texts: PyQGIS script works in QGIS console but produces no error or output as standalone

texts:
I'm trying map rendering with stand alone script, so I have small script which works fine in python console provided QGIS, however when I run same code in stand-alone script it executes without error, but output is not produced.

Here is my standalone script:

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtXml import *

QgsApplication.setPrefixPathddd("C:\Program Files (x86)\QGIS Valmiera", True)

QgsApplication.initQgis()

mapinstance = QgsMapLayerRegistry.instance()
mapinstance.removeAllMapLayers()

layer = QgsVectorLayer("F:\scripts\map\wb.shp", "test", "ogr")

if layer.isValid():
    print "shplayer succesfully loaded - adding to mapinstance"
    mapinstance.addMapLayer(layer)
else:
    print "shplayer failed to load!"

print mapinstance.mapLayers()

# create image
img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)
# set image background color
color = QColor(255,255,255)
img.fill(color.rgb())
# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)
render = QgsMapRenderer()
# set layer set
lst = [ layer.id() ] # add ID of every layer
render.setLayerSet(lst)
# set extent
rect = QgsRectangle(render.fullExtent())
#rect.scale(1.39800703)
rect.scale(1.1)
render.setExtent(rect)
# set output size
render.setOutputSize(img.size(), img.logicalDpiX())
# do the rendering
render.render(p)
p.end()

# save image
img.save('F:\scripts\out.png',"png")

QgsApplication.exitQgis()


texts: Using sf to combine polygons that share borders?

texts:
I have rows of respondents with their residential census tract polygon and am trying to generate first-order aggregated areas. Meaning I want a given respondent’s census tract polygon to be merged with all the polygons touching its borders. A separate file contains all possible US census tract polygons.

st_union seems appropriate only for aggregating together ALL interesting polygons to a single mass polygon.

Using st_intersection and then aggregating by row ID seems viable, but I’m unsure how to do the aggregation procedure.

I hope to remain within the sf environment.

Similar threads I’ve quered:
Dissolve only overlapping polygons in R

Joining polygons in R


texts: Loading Multiple CSV Files into QGIS 2.10

texts: How can I load multiple CSV files into QGIS with Longitude/Latitude already included in it.  I downloaded the "Load Them All" plug-in, however, the result loaded was converted as text files and cannot be visualize on the map canvass.  Any help? Are there any available plug-ins that can load multiple CSV file that will work fine in QGIS?


texts: Exporting to PDF: Multiple feature classes with similar names as separate PDFs

texts:
I have a list of feature classes used as facilities in network analysis to find out walking and cycling time. So for example, schools, hospitals, pharmacies etc. I therefore want to export a PDF with the schools plus the study area, hospitals plus the study area etc. as separate PDFs for each facility and associated analysis:

import arcpy

#Get the MXD
m = arcpy.mapping.MapDocument("CURRENT")

#Get the data frame
d = arcpy.mapping.ListDataFrames(m, "")[0]

#Get the layer list in the data frame
ls = arcpy.mapping.ListLayers(m, "", d)

#Set layer dictionary
layer_dict = {"layer_name" : ["turn_on_lyr_1", "turn_on_lyr_2"]}

#Get the folder you want to save your PDF files into
PDFPath = r'C:\Users\Charlie-x.Churchill\Documents\PyTrial'

#Create empty list
layerList = []

#add layer names to the list
for l in ls:
    layerList.append(l.name)

for lyr in arcpy.mapping.ListLayers(m, '',d):
    for lname in layerList:
      if lname in layer_dict:
           for    lyr in layer_dict[layer_name]: lyr.visible = True
           lyr.visible = True
           arcpy.RefreshActiveView()
           arcpy.mapping.ExportToPDF(m, PDFPath + "\\" + lyr.name + ".PDF")
           lyr.visible = False

del m


texts: How to convert 0 values to NaN values with rasterio

texts: I clipped my tif image to smaller tif via using a shapefile as mask with rasterio python library. After that my square shaped clip has bigger extent than shapefile and outside the the shp area is black and valued as 0. How can I change it to NODATA values?


texts:
I have been working with some cartodb infowindows and some dynamic templating in showing the infowindows when hovering. Sadly when hovering on different sublayers, interactivity of the different sublayers seems to trigger non related events on other sublayers or not trigger the infowindow at all. Any pointers on what is wrong here?

http://bl.ocks.org/jmwenda/17c06e7257963b3d0cee9a4ce84882ac


texts: Increasing label size of World_Street_Map layer or any other basemap layer in ArcMap?

texts:
I want to increase the inbuilt txt size of a world street basemap layer.

How to increase a txt size of a basemap layer in ArcMap?


texts: Renaming shapefiles using loop in Python

texts:
I'm trying to clean up a list of shapefiles that I've created using GME. I have a series of shapefiles (polylines and polygons) that are named "ID_final_kde.img_LSCV.shp" and "ID_final_kde.img_LSCV_poly.shp", where ID is unique to each individual. I'm trying to rename the polylines to "ID_final_LSCV.shp" and the polygons to "ID_final_LSCV_poly.shp", to no avail.

Here's what I've been trying:

# import modules
import arcpy, os, arcinfo, xlwt,time, glob
from arcpy.sa import *

arcpy.CheckOutExtension('spatial')

path = os.getcwd()
arcpy.env.workspace = path + '\\fawn_locations\\FINAL'
arcpy.env.overwriteOutput = True

lines=arcpy.ListFeatureClasses(feature_type='Polyline')
for each in lines:
    arcpy.Rename_management(each, each[:-16] + 'LSCV')

poly=arcpy.ListFeatureClasses(feature_type='Polygon')
for each in poly:
    arcpy.Rename_management(each, each[:-20] + 'LSCV_poly')

This code renames the .shp file, but not the associated .cfg, .prj, or .shx files, so I can no longer open the shapefile in ArcMap. How can I loop through the folder, renaming all the files while preserving their extensions?

I'm just learning python.


texts: In a versioned geodatabase, what impact do delta tables and the state tree have on query performance?

texts:
We have a versioned arcsde geodatabase (arcgis 9.3.1 on oracle 10g) with a fairly complex data model that includes about 100 featureclasses and non spatial tables, a geometric network and many relationship classes.

The data is edited daily by 5 or 6 arcmap users utilizing sde versioning. In addition versions are created by automatic services that interface with other business systems to perform edits in the geodatabase.
Query performance degenerates noticeably during the course of the day, so we have implemented a nightly script to achieve a full compress. On occasions when a relatively large number of edits are performed, the system can become unusable until after a full compress.

Its been suggested that oracle as configured can not come up with decent execution plans when confronted with these volatile delta tables. Is this a reasonable explanation? What approach should be taken to resolve it?

Update in response to comments

By end of day, the state tree is very linear, with only a little branching.
We compress nightly (get a full compress by deleting all versions).
Business tables are analyzed regularly.
Delta tables are not analyzed. They are locked (Attempt to analyze returns error "ORA-20005 object statistics are locked"). Neither are the volatile tables in the sde
schema - STATES, STATE_LINEAGES.


texts: What is required to get OGR to work with SQL Server?

texts:
All,

I've been having heaps of trouble getting OGR to connect to our SQL Server development environment. (read about it here and here)

SQL Versions: 2008, 2012
OS's: Windows 7, Ubuntu 12.04

My questions are:

1) What permissions does the user need in order to establish a connection to a SQL Server database with OGR?

2) Does anything need to be loaded on the SQL Server machine or application server itself?

I used to be able to connect even with QGIS, but that isn't working now either. (or here either)

Also: if you have a password with a ! in it, it will throw back a 'bash' error...fyi!

Ok thanks again, all!


texts: Why exporting map with qgis2web to Leaflet ends with errors?

texts:
Exporting qgis2web map with OpenLayers completes with no error.. but with Leaflet it ends with following errors..

An error has occurred while executing Python code:

IndexError: 1
Traceback (most recent call last):
  File "/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/maindialog.py", line 350, in saveMap
    feedback=self.feedback)
  File "/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/leafletWriter.py", line 101, in write
    folder=dest_folder)
  File "/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/leafletWriter.py", line 272, in writeLeaflet
    useVT, useShapes, useOSMB)
  File "/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/leafletLayerScripts.py", line 127, in writeVectorLayer
    useShapes, feedback)
  File "/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/leafletStyleScripts.py", line 65, in getLayerStyle
    sl, useMapUnits, feedback)
  File "/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/leafletStyleScripts.py", line 180, in getSymbolAsStyle
    sl = symbol.symbolLayer(sl)
IndexError: 1

Python version: 3.7.3 (default, Nov 10 2011, 15:00:00) [GCC 8.3.0]
QGIS version: 3.10.3-A Coruña A Coruña, 0e1f8464

Python Path:
/app/share/qgis/python
/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python
/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins
/app/share/qgis/python/plugins
/usr/lib/python37.zip
/usr/lib/python3.7
/usr/lib/python3.7/lib-dynload
/usr/lib/python3.7/site-packages
/app/lib/python3.7/site-packages
/app/lib/python3.7/site-packages/numpy-1.14.3-py3.7-linux-x86_64.egg
/app/lib/python3.7/site-packages/MarkupSafe-1.0-py3.7-linux-x86_64.egg
/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python
/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/qgis2web
/home/user/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins
/home/user/Documents/Web Map

Why it is happening??


texts:
Is there a way to sync GeoJSON files to my cartoDB table?

If this is not yet supported, what would be the best way of implementing this feature (syncing the data frequently with the existing cartoDB table) myself?


texts:
I use SRID 4326(WGS 84).
So what is the measure of distance in ST_ClusterWithin . m, km, degree ?
I ask it because in documentation it only says:

geometry[] ST_ClusterWithin(geometry set g, float8 distance);

while ST_DWithin function signature at the same time looks like:

boolean ST_DWithin(geometry g1, geometry g2, double precision distance_of_srid);

boolean ST_DWithin(geography gg1, geography gg2, double precision distance_meters);

boolean ST_DWithin(geography gg1, geography gg2, double precision distance_meters, boolean use_spheroid);

As i know distance_of_srid for 4326 is degrees.


texts: ArcCatelog 9.3 feature classes with annotation relationships allowed in a geometric network

texts: I have a gas utility dataset that we are working on getting converted into a geometric network for easier analysis.  The feature classes are part of an ArcSDE and have a ton of relationship classes attached, mostly annotations. I know that geometric networks will not accept feature to feature relationship classes when building the network, but what about relationship classes built from annotations? If we cannot get those annotation classes into the geometric network, then we might end up scrapping the idea all together in our department.


texts:
I am trying to join points to polygons. To be more precise I have a dataset with points associated with a particular political event. Each point comes with latitude and longitude and other information. And I have a shapefile that has polygons that have regional data information at the Nuts1 level.

I want each of the points to be given the NUTS code in which it belongs to based on latitude and longitude.

I do spatial join on ArcGIS joining points to polygons. I have the same projection in both files. But when the join is completed the data from the polygons do not appear in the joint file. I have the corresponding columns, e.g., Nuts code, or population of this NUTS but all info comes with 0's.

Any suggestions why this happens and how can I fix it?


texts: Re-using dataframe in ArcGIS without ArcMap crashing with R6025 error?

texts:
I have a large map with an inset. In Arc9.3 I was able to copy that inset data frame into another mxd to use again as another inset - it is set up exactly how it is needed. Rather than re-build a new dataframe and add all the layers and symbolize, its way easier to re-use a dataframe already created. When I try to do this in Arc10.0, the dataframe copies okay however when I try and access the dataframe properties, it crashes ArcMap. (The only adjustment I need to do in the properties is to reset the map extent frame.) The steps I took to access properties, right click & choose properties. Upon chosing properties, ArcMap crashes. This process worked so smoothly and without error in previous versions. I feel frustrated with this simple task. Is there another way to re-use an existing dataframe, without ArcMap crashing?

Not every time, but sometimes this error pops up before the crash:


texts: Postgis - Zonal statistics - Parallel computation?

texts:
I am looking to speed up a postgis query, regarding zonal statistics. I am using postgis 3.0 with postgreSQL 12 and I was hoping that I could apply the below script to speed up processing using multiple CPUs;

// Parallel processing has no effect
SET max_parallel_workers = 8;
SET max_parallel_workers_per_gather = 4;

// Zonal statistics
DROP TABLE IF EXISTS schema.zonal_stats;
CREATE TABLE schema.zonal_stats AS
SELECT id,geom,
(ST_SummaryStats(ST_Union(ST_Clip(raster,geom)))).sum
FROM schema.raster, schema.ply
GROUP BY id,geom;

I could use the rasterstats and multiprocessing packages in Python, which does provide a computation speedup. However, I would like to apply parallel processing natively in postgis, as my data is stored in a postgreSQL db.


texts:
I have a problem with a shapefile I created exporting features from another shapefile.

I have all my lines visible but their rows only appear in the attribute table window when they're selected, do you know why and how can I fix this ? (I posted pictures to illustrate my problem)

selected lines

unselected lines


texts:
I edit my raster catalogs to get a specific symbology and then I have to manually export the XML. I then load the XML into the mosaic dataset. I can't seem to find a way to automate the process of exporting all the XMLs and then importing them into the mosaics. Has anyone ran into this issue before?

Also, I've looking into exporting the symbology to a layer and then pulling the stats from that but it failed. I almost thought the Create Referenced Mosaic Dataset tool was going to be exactly what I was looking for but that also was a dud.


texts: Show only current atlas feature (line) within current atlas plot

texts: I am using lines as atlas features so that I can rotate the view by the azimuth of the line.  How can I have only the current atlas feature printed on the current sheet?  I can see how "Within" would work for polygons but can't find a similar solution for a line feature.


texts: Classifcation Time Series Error in Google Earth Engine?

texts:
In the below code I was going to do classification process for time series data but in training data collection returned this error:

code link: https://code.earthengine.google.com/b81fe0a02ece289b3c788a5ccba67ac9
Map.centerObject(table);
Map.addLayer(table);

var start = '2018-01-01';
var end = '2019-01-01';

var persiannData = function(img){
  return img.clip(table)
  .multiply(255 / 2.0).toInt()
  .copyProperties(img,['system:time_start','system:time_end']);
};

var persiann25km = ee.ImageCollection("NOAA/PERSIANN-CDR")
.filterDate(start, end)
.filterBounds(table).map(persiannData);

print('persiann25km',persiann25km);

/// ancillary data

var dem = ee.Image("USGS/GTOPO30")
.clip(table);

var ndvi = function(img){
  return img.clip(table).addBands(dem)
  .copyProperties(img,['system:time_start','system:time_end']);
};

var evi = function(img){
  return img.clip(table)
  .copyProperties(img,['system:time_start','system:time_end']);
};

var lst = function(img){
  return img.clip(table).select('LST_Day_1km').multiply(0.02)
  .copyProperties(img,['system:time_start','system:time_end']);
};

var modisNDVI = ee.ImageCollection("MODIS/MOD09GA_006_NDVI")
.filterDate(start, end)
.filterBounds(table)
.map(ndvi);

var modisEVI = ee.ImageCollection("MODIS/MOD09GA_006_EVI")
.filterDate(start, end)
.filterBounds(table)
.map(evi);

var modisLST = ee.ImageCollection("MODIS/006/MOD11A1")
.filterDate(start, end)
.filterBounds(table)
.map(lst);

// data integration

var modisDataset = modisNDVI.combine(modisEVI).combine(modisLST);

print('modisDataset',modisDataset);

var innerJoin = ee.Join.inner();

var filterTimeEq = ee.Filter.equals({
  leftField:'system:time_start' ,
  rightField: 'system:time_start'
});

var innerJoinModis = innerJoin.apply(modisDataset, persiann25km, filterTimeEq);

print('persiann & modis', innerJoinModis);

var datasetMap = ee.ImageCollection(innerJoinModis.map(function(feature){

  return ee.Image.cat(feature.get('primary'), feature.get('secondary'));

}));

print('datasetMap',datasetMap);

// modelling process

var classificationProcess = datasetMap.map(function(img){

  var bandNames = img.select('NDVI','EVI','elevation','LST_Day_1km')
  .bandNames();

  var trainingData = img.stratifiedSample({
    numPoints: 100,
    classBand: 'precipitation',
    region: table,
    scale: 1000,
    });

  var classifier = ee.Classifier.smileRandomForest(80)
  .train({
    features: trainingData,
    classProperty: 'precipitation',
    inputProperties: bandNames
  });

  var classified = img.select('NDVI','EVI','elevation','LST_Day_1km')
  .classify(classifier);

  return classified
  .copyProperties(img,['system:time_start','system:time_end']);

});

print('classificationProcess',classificationProcess);


texts: Load Vector Layer from table using processing in QGIS3

texts:
I am trying to turn a CSV of data into a shapefile using processing.

mStasIn = "C:\GIS_Abschluss\AirData\MeasurementStations.csv"
mStasOut = "C:\GIS_Abschluss\AirData\MeasurementStations.shp"

mStas = QgsVectorLayer(mStasPath, "mStas", "ogr")

processing.run('qgis:createpointslayerfromtable',
            {'INPUT' : mStasIn,
            'XFIELD': 'xcoord',
            'YFIELD': 'ycoord',
            'ZFIELD': 0,
            'MFIELD': 0,
            'TARGET_CRS': 'EPSG:32633',
            'OUTPUT': mStasOut})

iface.addVectorLayer(mStasOut, 'mStas', "ogr")

This creates a shapefile, but there are no points on my canvas.


texts: Change maxZoom level on Leaflet Marker Cluster Group?

texts:
Our leaflet web app allows users to switch between tile layers.  The tile layers sometimes have different max zoom implemented by the tile provider or local config.

We have another layer with markers, clustered using markerClusterGroup configured to open up the clustering at max zoom AKA spiderfy.  The problem is that when the user switches to a tile layer with a smaller max zoom the spiderfy function no longer works.

new L.markerClusterGroup({ spiderfyOnMaxZoom: true }) ...

How can we dynamically tell the markerClusterGroup that max zoom level changed ?


texts: Can draw points but not polylines in WPF Runtime code

texts:
The following C# code works when drawing in Point mode, but not Polyline . . .

  // "Draw" is ESRI.ArcGIS.Client.Draw object
       _draw = new Draw(_map);
       _draw.DrawMode = DrawMode.Polyline;
       _draw.DrawComplete += (senderObject, drawEventArgs) =>
        {
            _barrierLine = new Graphic();
            _barrierLine.Symbol = this.Resources["LineSymbol"] as SimpleLineSymbol;
            _barrierLine.Geometry = drawEventArgs.Geometry;
            _lineGraphicsLayer.Graphics.Add(_barrierLine);
        };

This is what the XML code looks like (I tried several different methods):

   <Window.Resources>
        <esri:SimpleMarkerSymbol x:Key="PointSymbol" Color="DarkViolet" Size="15" Style="Diamond"></esri:SimpleMarkerSymbol>
        <esri:SimpleLineSymbol x:Key="LineSymbol" Color="Aqua" Width="5" Style="Solid"></esri:SimpleLineSymbol>
    </Window.Resources>
    <Grid>
        <Grid.Resources>
            <esri:SimpleMarkerSymbol x:Key="RedMarkerSymbol" Color="Red"/>
            <esri:SpatialReference x:Key="theSpatialReference" WKID="4326"/>
            <esri:SimpleLineSymbol x:Key="LineSymbol" Color="Aqua" Width="3" Style="Solid"/>
        </Grid.Resources>
            <!--Graphics Layer-->
            <esri:GraphicsLayer ID="graphicsLayer" x:Name="_graphicsLayer">
                <esri:Graphic>
                    <esri:Graphic.Symbol>
                        <esri:SimpleMarkerSymbol Style="Circle" Color="Red"/>
                    </esri:Graphic.Symbol>
                    <esri:Graphic.Geometry>
                        <esri:MapPoint  X="0" Y="0" >
                            <esri:MapPoint.SpatialReference>
                                <esri:SpatialReference WKID="4326"/>
                            </esri:MapPoint.SpatialReference>
                        </esri:MapPoint>
                        </esri:Graphic.Geometry>
                </esri:Graphic>
            </esri:GraphicsLayer>
            <esri:GraphicsLayer ID="lineGraphicsLayer" x:Name="_lineGraphicsLayer">
                <esri:Graphic>
                    <esri:Graphic.Symbol>
                        <esri:SimpleLineSymbol Style="Solid" Color="Aqua" Width="3" />
                    </esri:Graphic.Symbol>
                    <esri:Graphic.Geometry>
                        <esri:Polyline>
                            <esri:Polyline.SpatialReference>
                                <esri:SpatialReference WKID="4326"/>
                            </esri:Polyline.SpatialReference>
                        </esri:Polyline>
                    </esri:Graphic.Geometry>
                </esri:Graphic>

            </esri:GraphicsLayer>
        </esri:Map>
        <StackPanel Orientation="Vertical" HorizontalAlignment="Left">
            <Button Content="Stops" VerticalAlignment="Top" Name="StopsButton" Margin="5,5,5,0" Click="StopsButton_Click" />
            <Button Content="Point Barriers" Height="23" Name="PointBarriersButton" Width="76" Click="PointBarriersButton_Click" />
            <Button Content="Solve" Name="SolveButton" Margin="5,5,5,0" Click="SolveButton_Click" />
            <Button Content="Line Barriers" Height="23" Name="LineBarriersButton" Width="76" Click="LineBarriersButton_Click" />
        </StackPanel>
    </Grid>


texts: Spatial join string attributes for intersecting polygons (QGIS)

texts:
I have seen tutorials for the 'Spatial Join' tool that allows me to summarize statistical data for overlapping polygons on two separate layers. But I simply want to append a 'String' attribute from an underlying layer to polygon features on an overlying layer.

Specifically: I have a layer of 'land uses'; and a layer of 'buildings'. I would like to attach the 'land use [string]' attribute for a particular 'land use' polygon (on the land-use layer) to all the buildings (smaller polygons on the 'buildings layer') that intersect with that specific land use polygon.

Can somehow tell me what tool or script to use with QGIS? (And how to use it.)


texts: Create grid from raster with equal size squares

texts:
I have a raster of the globe with coordinate system WGS_84.

I want to create a grid on that raster where every square has the same size.

When I use the Create Fishnet tool the size of squares get bigger as I move away from the equator.

Is there a way to make a grid where every square has the same size or at least where the differences between squares are small?


texts: OSMDroid app using custom ArcGIS Tile Server - tiles are shuffled

texts:
I am developing an android mapping app using OSMDroid.  I am attempting to use free custom aerial imagery, completely independent of Google and/or Bing APIs.  Please, do not propose any solution that uses their mapping APIs.

I have managed to display satellite imagery by including this code:

mapView.setTileSource(TileSourceFactory.MAPQUESTAERIAL);

but Tile Server does not offer tiling above 11 zoom and I need to get closer than that (say 15-16?).

Using ArcGIS tile server, I manage to display satellite imagery even to 16 layer zoom level, but tiles are shuffled around.

mapControl = (MapController) mapView.getController();
mapControl.setZoom(11);
String[] urlArray = {"http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/"};
mapView.setTileSource(new XYTileSource("ArcGisOnline", null, 0, 18, 256, ".png",urlArray ));

Basemap tiles are shuffled and do not correspond to lat/lon, but overlay is ok.


texts: Differentiating between children and adults depending on AGE using Cartodb Torque wizard?

texts:
I'm trying to differentiate between children and adults depending on AGE, The age of these individuals are contained under the Column "palestinian_civilians" .I'm using Carto Torque wizard.Below is a little code I was using.Nothing is happening, what am I doing wrong?

torqueLayer.setSQL("select *,(CASE WHEN palestinian_civilians ='0' THEN 1 ELSE 2 END) as palestinian_civilians_n from table");

Map {
-torque-frame-count:512;
-torque-animation-duration:30;
-torque-time-attribute:"palestinian_civilians";
-torque-aggregation-function: "round(avg(palestinian_civilians))";
-torque-resolution:2;
-torque-data-aggregation:linear;
}
#layer {
...
  marker-fill: #FF0000;
  // avg of 0 and 36
  [value > 18] { marker-fill: #0000FF; }
...
}


texts: WMS of Geoserver SQL view shows nothing, however WFS works properly

texts:
I am using Geoserver to create a WMS service to serve data out of an SQL Server database. Everything works fine (WMS & WFS) untill I create an SQLView in Geoserver.

The SQLView:

SELECT *
FROM rmaster

I have defined the CRS and geometry column. The WFS service of the SQLView works properly and serves the data, however the WMS serves empty images. Does anyone have any idea what's causing this?


texts:
I wondered if anyone could help. I have setup osrm on a spare computer it is running on Ubuntu 14. It works fine and I can query it using the default car setting for transport. I would like to be able to setup and use different modes of transport ie foot and public transport as well. the query I send to osrm is this:

$json = file_get_contents("http://192.168.0.4:5000/viaroute?loc=$startlat,$startlong&loc=$endlat,$endlong&instructions=false&alt=false&geometry=false");


texts: I am trying to figure out a way to create a stream order file with the Horton method, however seeing as Arcmap only does either Strahler or Shreve, this isn't ideal. I opted to try QGIS, however it cannot apply the Fill sinks algorithm on the World 30 DEM over the state of catarina (or even a smaller sized region). I've disabled all unnecessary plugins and only have "Processing" enabled. It constantly freezes. I figured it was simply loading however after a night of "working", it has failed to display the information.


texts: Convert exported attribute table back to shapefile in QGIS

texts:
I exported the attribute table of a shapefile to Excel (also saved in CSV). And I'm trying to add some columns to the table, then convert the table back to shapefile with the new columns. The shapefile is a network of roads, but I wasn't able to find a columns that contains x/y coordinates or WKT in the attribute table of the original shapefile.
I'm wondering how can I convert this exported and modified attribute table back to shapefile? I don't have ArcGIS so I think I have to use QGIS or other method.


texts: Put color for raster data(aspect.tif) on QGIS1.8.0

texts:
I create aspect.tif from DEM file.
I could put the values(0 to 360)from contrast enhancement.

I could set color from color map(freakout or Pseudocolor).
In here, the data is shown by only 4 colors.
I could not change colors and add color type.

First of all, I could not select colormap tab.

My ideal color setting is 8 colors with 8 direction(0 to 45, 45 to 90, ....315 to 360)

Please tell me how to change color.


texts: Batch Importing 3.0ft Esri Grid to QGIS?

texts:
How to batch import Esri Grid files to QGIS?

I have received Bare Earth DEM files for a large swath of land. Each file has a unique identifier referencing a certain ESRI grid (eg. be_ef_001,be_EF_002, etc). Each grid is 1/2 mile X 1/2 mile.

In order to import each grid in QGIS I open the selected folder and select the .adf file. There are over 1,000 files that I am working with and it seems that I have to do this individually for each one.

ArcGIS Desktop can view each folder as a directory and I can simply select all and import. I do not have access to ArcMap or ArcGIS Pro, only QGIS. This is something that will likely be a reoccurring issues, and I would like to find a way where I don't have to import each grid manually.


texts: Zooming to extent of query result in ArcGIS API for JavaScript?

texts:
How do I change the extent on the map in function of a query?

in legend I  have no problem using the function showResults (featureSet) {...
but,  can not get any display in AMD. Here I put a piece of script that I want to switch from legend to AMD:

function showResults(featureSet) {
  //remove all graphics on the maps graphics layer
  map.graphics.clear();
  //QueryTask returns a featureSet.  Loop through features in the featureSet and add them to the map.

  //Performance enhancer - assign featureSet array to a single variable.
  var resultFeatures = featureSet.features;
  var extent;
  for (var i = 0, il = resultFeatures.length; i < il; i++) {
    //Get the current feature from the featureSet.
    //Feature is a graphic
    var graphic = resultFeatures[i];
    graphic.setSymbol(symbol);
    geometry = graphic.geometry;
    ext = geometry.getExtent();
    if (extent) {
      extent = extent.union(ext);
    } else {
      extent = new esri.geometry.Extent(ext);
    }

    //Set the infoTemplate.
    graphic.setInfoTemplate(infoTemplate);

    //Add graphic to the map graphics layer.
    //map.addLayer(stateLayer);
    map.graphics.add(graphic);
  }

  map.setExtent(extent, true);

}


texts:
I'm new to maps, and have been searching for the past couple of hours how to get started with OpenStreetMaps without much success.

What I would like is a spreadsheet with the coordinates, names and addresses of all the fast food restaurants in Dublin, Ireland which I could then paste into something like batchgeo and explore.

I initially tried to go to openstreetmaps and restrict what it showed to just the Tag:amenity=fast_food but I couldn't see a way to do that.

So I've now downloaded the osm file from cloudmade and opened it in JOSM, but I'm not seeing how to remove the layers.

Could one of you guys tell me the simplest way to achieve what I'd like? Thanks.

An update for those keeping score, I've now got the file into osmosis and am trying to filter it that way.


texts: Polygons and lines not drawing correctly (in QGIS from PostGIS?

texts:
I load osm file into postgresql and then i connect to the database with qgis t to plot it and  i see this

it happens to whatever file i import to the database, even when i download it from http://metro.teczno.com/
i also get the same result when i plot the lines but i did not include them in the screenshot just to be clearer [ i can provide one if asked ]...points are fine...

some relevant info that might be helpful

xubuntu 12.04
postgresql-9.3-postgis-2.1

steps:

osm=#create extension postgis;
osm=#\i /usr/share/postgresql/9.3/contrib/postgis-2.1/legacy.sql
osm=#\i /usr/share/postgresql/9.3/contrib/postgis-2.1/legacy_gist.sql

then i simply load the osm file with

postgres@desk: osm2pgsql -U postgres  -d  osm  -S /usr/share/osm2pgsql/default.style mapFile.osm

output seems fine (at least to me)

Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
Setting up table: planet_osm_line
NOTICE:  table "planet_osm_line_tmp" does not exist, skipping
Setting up table: planet_osm_polygon
NOTICE:  table "planet_osm_polygon_tmp" does not exist, skipping
Setting up table: planet_osm_roads
NOTICE:  table "planet_osm_roads_tmp" does not exist, skipping
Mid: Ram, scale=100

!! You are running this on 32bit system, so at most
!! 3GB of RAM can be used. If you encounter unexpected
!! exceptions during import, you should try running in slim
!! mode using parameter -s.

Reading in file: mapFile.osm
Processing: Node(13k) Way(2k) Relation(258)  parse time: 1s

Node stats: total(13746), max(2147483647)
Way stats: total(2090), max(267296233)
Relation stats: total(258), max(3555834)

Writing way (2k)

Writing relation (258)
Committing transaction for planet_osm_line
Committing transaction for planet_osm_point
Committing transaction for planet_osm_roads
Sorting data and creating indexes for planet_osm_line
Committing transaction for planet_osm_polygon
Sorting data and creating indexes for planet_osm_point
Sorting data and creating indexes for planet_osm_roads
Sorting data and creating indexes for planet_osm_polygon
Completed planet_osm_roads
Completed planet_osm_line
Completed planet_osm_polygon
Completed planet_osm_point

Can someone explain why this might happen?

As far as i can see, some polygons, if not all, polygons are not correctly parsed and closed and include a point that is located at the center of the plot(?)

I am really new to this and my first goal is just to correctly load an osm file into a postgis database and just display it correctly (i actually want to render it with python mapnik, but that will follow if/when this is solved).


texts: Splitting multiline into several same-length multilines

texts:
Is there a method to convert a multiline geometry into several multiline geometries having the same length?

For example break this:

MULTILINESTRING((1 1, 3 1, 3 1.5, 2.5 1.5))

Into [approximately] same-length segments with the length of 1:

MULTILINESTRING((1 1, 2 1))
MULTILINESTRING((2 1, 3 1))
MULTILINESTRING((3 1, 3 1.5, 2.5 1.5))


texts: Network Analyst for QGis or another Open Source Geographic Information System?

texts:
Im looking for alternative ways for analyse my network dataset.
So, there is another way to analyse an GIS network? I mean if there are any useful plug-in capable to determine service areas, routes, or even VRPs like Arcmap does? Road Graph plugin can do all this stuff?


texts: There is a way to merge shp polygon with different fields? or there is a tool to copy the fields and format from one shp to another (not the content)?


texts: Updating MapGuide Enterprise 2011's Query Tool

texts:
I am trying to create a custom theme for my company's website using MapGuide Enterprise 2011, and I am looking for advice on fine-tuning the software's built-in query tool to allow for more advanced querying.

The built-in Query tool is (as far as I can tell) only able to perform a single query on a layer, however I am hoping to be able to query multiple layers/multiple queries in the same layer for more precise searching. Being able to save my query results as a selection would also be ideal.

Any assistance/feedback would be greatly appreciated. Thank you.


texts: Georeferencing US County names in CartoDB

texts: I have a massive file with county-level data from the 19th century. I'd like to georeference it, but I'm not sure what I need to do, as CartoDB isn't recognizing the county names through the 'Admin. Regions.' Do I need to import shapefiles and somehow connect them?


texts: Calculation GPS Accuracy based on known point

texts:
I have 3 GPS's that were on the same point for about an hour.
Each of them yield a fairly accurate position but not a complete accurate location (obviously)

I want to know how to calculate the best GPS out of the three. The information I have is the X Y Z of the actual point, and for each GPS I have (TIME , X, Y, ACCUARCY)

What I did so far is:

basic distance error calculation - Average of each gps and compare
standard deviation of errors and knowlocation - compare


texts: Using same ROI/Signature for classifying spatially different rasters?

texts:
Using QGIS SCP Plugin, I have classified a raster of a state.

Can I use the same ROI/signatures to classify rasters of other states?


texts: Is there a significantly faster way to buffer than arcmap?

texts:
I hope this question is appropriate here and I apologize if it is more of a discussion question than a specific technical problem.

I have to buffer 150,000 polygons for a real estate project.  While it's not impossibly slow, it got me wondering about 'the fastest way to geoprocess'.
What packages, tips, or tricks could be used to make this faster?
Does anyone have an idea how arcmap, arcpy, shapely, postgres compare in performance time?

EDIT: After asking around, I heard 'arcpy is slow!' a lot.


texts: Removing no data area (sea) from raster in ArcMap?

texts:
I have some difficulties removing the no data part from a GeoTIFF SRTM. I need to remove all sea parts. I'm new in ArcGIS.

This is the result i wanted


texts: How to prepare SVG images for presentation in OpenLayers?

texts: I have a collection of SVG vector images I need to display on a map, probably using OpenLayers/KML. Is there any way to convert a vector image from SVG to KML? Or would it make more sense just to use PNG?  I would prefer to use vectors because the images will be rotated and need to scale properly when zooming.


texts: Stacked features on same layer, can't bring on foreground

texts:
So I have multiple features on the exact same coordinates (which happen to be 0,0 but it is irrelevant I believe)

I run a routine that does something like:

clientVectorSource.forEachFeature(function(item) {

                        if (item.get('id') == 1233){
                            /* I tried this part
                            var tempItem=item;
                            clientVectorSource.removeFeature(item);
                            clientVectorSource.addFeature(tempItem);
                            */
                            map.getView().setCenter(item.getGeometry().getCoordinates());
                            selectedFeature=item;
                            if (map.getView().getZoom() < 15)
                                map.getView().setZoom(15);
                        }
                    });

selectedFeature is a global variable not relevant to my problem.

The problem is that the "item" is sometimes stacked below the rest of the features due to how it happened to be rendered at the time of execution. I want to bring it on the foreground, but I can't.

I know it could possibly be solved by changing layers on the feature and then changing zIndex to the entire layer but I don't want to resort to that since it is not a very elegant solution and can possibly lead to problems in the future.

I also tried changing the style of "item" to the following and it didn't work:

var iconStyleNew = new ol.style.Style({
                    image: new ol.style.Icon( /** @type {olx.style.IconOptions} */ ({
                        anchor: [0.5, 30],

                        anchorXUnits: 'fraction',
                        anchorYUnits: 'pixels',
                        src: 'data/iconMarker.png'
                    })),
                    zIndex : 999

                });


texts:
I'm trying to split polygons into equal areas. While fishnet sort of works; the areas I'm trying to split don't touch; they're national parks. Fishnet spans the whole area, not just the polygons  in question. Ultimately, I'm searching for a method to split polygons into a user-given number of equal area polygons.

I've read everything I can about the Huber methods to do this in ArcView, but without access to that (or experience with it) this isn't the answer for me.

Thanks very much.