texts: Additional Basemap Menu in ArcGIS Viewer for Flex?

texts:
I am using ArcGIS Viewer for Flex 3.5 to create an application. As it covers more than one areas of interest, I need to provide an option to the user to only select base maps from his own desired domain. In order to do that,  I need more than one Basemap dropdown menus in the application whereas, the default map switcher widget only provides options for one menu.

I have not been able to find any worthwhile support or guidance for the same. The only relevant source of information which I could find is a thread related to similar task in sample flex viewer in old ESRI forums at this link:
http://forums.esri.com/Thread.asp?c=158&f=2421&t=282643.

However, the flex code when compared the 3.5 version is no more relevant and I cannot figure out how to make use of that.


texts: How to find peaks by comparing with neigbhour points

texts:
I have a raster layer with 3 bands
band1 - elevation data
band2 - x-coord
band3 - y-coord

I want to run  a grid of 9*9 cells over this raster layer staring from top-left to bottom-right row-wise and to get the peak (highest altitude) out of the neighboring cells in the current 9*9 grid and put that peak into a 2D-matrix which contains XY coordinate of the corresponding peak after completing the iteration, I will be left with a matrix of coordinates I came to know that grid function from raster package in R will be useful any more methods.


texts: Calculating grid convergence (True North to Grid North)

texts:
I would like to know how to calculate the grid convergence angle between True North and Grid North for UTM maps.
I am using ArcGIS 10.0 (ArcView licence). and have used the Calculate Grid Convergence Angle tool from the Cartography toolbox.
However, I would like to know how to calculate this angle myself (using Excel) to compare against the ArcGIS result?


texts: Using Multiple Raster Intermediate Inputs in Raster Calculator within QGIS Process Modeler

texts:
I am stringing together a series of raster calculations in the QGIS process modeler in order to calculate a surface temperature dataset from Landsat and weather data.

When I use the output from multiple prior algorithms in the model as the input for a raster calculation, the result is a blank grid (-1.79769e+308 to 1.79769e+308). Calculations performed with input parameter bands work as expected.

This issues appears to have been addressed in this QGIS bug report.

However I've updated my install and the problem persists. Can anyone advise on this? I'd rather not split my model into several smaller models.


texts: Rotating orthophoto based on polygon boundary (rotated pixel)?

texts:
When I clip an orthophoto using a vector polygon, the resulting pixels orientation is north-up like as shown in the graphic below. I used Global Mapper, ERDAS, Arcmap "Clip" "Extract by Mask", "Extract by polygon" and all of them outputs North-up pixels which ends up in big files

What I am trying to accomplish is to rotate the pixels based on the polygon boundary as shown in the image below.This task is accomplished by software called SOCET SET and the result is a much smaller file in size.

Does anybody know if there is another software besides SOCET SET that can rotate the pixel based on a rotated rectangular vector polygon?


texts: Matching and aggregating data between two overlapped layers using QGIS

texts:
I'd like to know how a vector layer can take and aggregate attributes of a field, of another vector layer which is overlapped.
For example, a layer T1 has not overlapped records while on the contrary the layer T2 could have overlapped records. I should have a "NEW" field in T1 with the data of T2, separated with ';' every time T2 records overlaps on T1.

How do I automate this process using QGIS?


texts: How to find Drawn Polygon intersection with layer elements in OpenLayers

texts:
I draw polygon vectors and I can find intersection to each other:

vactorTest.features[0].geometry.intersects(vactorTest.features[1].geometry)

but I want to find which layer elements this polygon intersects.

How can I find it?

EDIT:

I have vector layer like this:

layerArray = new OpenLayers.Layer.WMS(wms_layer_configs[i].label, "http://somelink/geoserver/wms",
        {
        // srs: 'EPSG:xxxx',
            layers: workspace + ":" + wms_layer_configs[i].name,
            tilesOrigin: map.maxExtent.left + ',' + map.maxExtent.bottom,
            format: 'image/png',
            transparent: true
        }, {
            opacity: 0.75,
            isBaseLayer: false,
            visibility: false
        });

than I draw polygon and want to display on the map features from layer which intersects this drawn polygon


texts: How to integrate XML markers with MarkerCluster?

texts:
I'm having a problem integrating XML markers with MarkerClusterer.
I need to integrate these two and show only the markers within bounds. I'm able to load the results in the XML to the map, but can not integrate the MarkerClusterer.

var customIcons = {
      Commercial: {
        icon: 'http://cdn1.iconfinder.com/data/icons/free-business-desktop-icons/32/Company.png',
      },
      Condo: {
        icon: 'http://cdn1.iconfinder.com/data/icons/REALVISTA/project_managment/png/32/agency.png',
      },
      Home: {
        icon: 'http://cdn1.iconfinder.com/data/icons/free-business-desktop-icons/32/Home.png',
      },
      Lot: {
        icon: 'http://cdn1.iconfinder.com/data/icons/Real-Estate-png/32/Lands.png',
      }
    };

function load() {
  var map = new google.maps.Map(document.getElementById("map"), {
    center: new google.maps.LatLng(<?php echo $lat ?>, <?php echo $long ?>),
    zoom: <?php echo $zoom ?>,
    mapTypeId: 'roadmap'
  });
  var infoWindow = new google.maps.InfoWindow;

  downloadUrl("mapxml.php", function(data) {
    var xml = data.responseXML;
    var markers = xml.documentElement.getElementsByTagName("marker");
    for (var i = 0; i < markers.length; i++) {
      var name = markers[i].getAttribute("name");
      var address = markers[i].getAttribute("address");
      var type = markers[i].getAttribute("type");
      var point = new google.maps.LatLng(
          parseFloat(markers[i].getAttribute("lat")),
          parseFloat(markers[i].getAttribute("lng")));
      var html = "" + name + "
" + address;
      var icon = customIcons[type] || {};
      var marker = new google.maps.Marker({
        map: map,
        position: point,
        icon: icon.icon,
      });
      bindInfoWindow(marker, map, infoWindow, html);
    }
  });
}


texts: Is there a character I need to use to state a new line in the MapInfo SQL Select box (e.g like ; in SQL Management Studio) or is this not possible. The highlighted (yellow) area is where the new line begins.


texts: error: wrong resolution for raster calculator, how to work with this?

texts:
I have two raster files (lets call them raster A and B) and I would like to multiply them with the raster calculator. I keep getting an error message when working with the raster calculator. I thought it could be related to differences in resolution of the rasters. Raster A has a cell size of 25x25m and raster B of  0.00029444091x 0.00029444091 m.

I wanted to change the small cell size so I tried the aggregate tool and the resample tool. It did not work.
The resample tool gave the following error message:

Does someone has any ideas how to work with this? Are there any other tools available to change the resolution?

I am using ArcGIS 10.2.2


texts: What are problems & biases stemming from mathematics behind using Monte Carlo Markov Chain to make heatmap?

texts:
My buddy and I are using Monte Carlo Markov Chains to make a heat map - in the interest of rigour, and out of curiosity:

"What are the problems & biases that stem from the mathematics behind using Monte Carlo Markov Chain to make a heatmap, and what alternatives or corrective measures would you use to remedy this?"

I'm finding it difficult to find a comparison of different predictive algorithms ! What I want to know is what biases are induced by the mathematics behind this algorithm, and which simpler/more complex proceedures would one use to correct for these.


texts: Understanding Join Attributes by Location in QGIS?

texts:
I wonder if there is somewhere an explanation for the parameters that can be set when running the Algorithm "Join Attributes by Location" in QGIS 2.18.

I had a look here:

Documentation QGIS Testing | Join attributes by location
Documentation for QGIS 2.8 | fTools Plugin

but there is no explanation. Any hints?


texts: Renaming shapefiles and raster files with Python

texts:
I'm trying to add "_Kiwi" to the name of multiple shapefiles and raster files in a folder C:\data.  For the files with multiple extensions (.shp.xml for example), I can't get the files to rename correctly.

My python script (in c:\data\renamefiles.py) is:

import glob
import os
import sys

dirToProcess = sys.argv[1]
matchPattern = sys.argv[2]
replacePattern = sys.arg[3]

print "Processing " + dirToProcess

if (os.path.isdir(dirToProcess)):
    print dirToProcess + " found"
else:
    print dirToProcess + " not found"
    exit ()
print "continuing processing..."

def rename (dir, pattern, titlePattern):
    for pathAndFilename in glob.iglob(os.path.join(dir, pattern)):
        title, ext = os.path.splitext(os.path.basename(pathAndFilename))
        print "found: " + title + "renaming..."
        os.rename(pathAndFilename,
                           os.path.join(dir, titlePattern % title + ext))

The files with multiple extensions are not getting renamed correctly.  Here are two examples:

Example 1:   original files is: shapefile1.shp.xml

If I type....  python renamefiles.py c:\data *.shp.xml %s_KIWI, I end up with: shapefile1.shp_KIWI.xml

Example 2:  original file is: Raster1.sid.aux.xml

If I type python renamefiles.py c:\data *.sid.aux.xml %s_KIWI, I end up with: Raster1.sid.aux_KIWI.xml

How can I do this with only my current libraries/modules?


texts: Where are coordinates stored in ArcGIS Pro?

texts:
I am brand new to ArcGIS Pro. I have never used this software before and know nothing about it.
I have a shapefile that I have imported to ArcGIS and I can't find the coordinates in the table. It has all these green dots that have been plotted all over the map.
I can click each dot and manually copy and paste the DMS coordinate data for each dot. But that's obviously a waste of time.
When I go into the table of the imported data, it has all types of information on each dot, but there is no field with the same coordinate data as shows up when I click the dot.
Where do I find this coordinate data?
I am simply looking to export to a readable file, a list of all the things that have been plotted on the map and what their coordinates are.


texts: Defining temporal extent of INSPIRE data when it changes every day?

texts:
I am trying to create INSPIRE compliant metadata for a view service.
I am struggling with temporal extent however. The data concerns projects which started being recorded back in 2010, however the data available to the service will only show data concerning projects who are executed in the last year or the next 6 months. Of course, this is updated each day, so my temporal extent is moving.

How can I record this in the metadata?


texts:
I have downloaded a digital soil map of the world from: https://worldmap.harvard.edu/data/geonode:DSMW_RdY

It has lots of files in it with different endings (shp, sbx, xls etc.). I am using QGIS 3.4. When I add vector layer and add the shapefile, all it does is loads a bunch of polygons with no soil colour to differentiate them (they are all the same colour). When I try to load all the files in the zip file at once, the same thing happens.

Can anyone tell me what I am doing wrong?


texts:
When I try to save a specific vector shapefile (line) as another  vector shapefile line also (Backup Copy) in QGIS version 2.10 or 2.18.3 I get the following error message:

  Export to vector file failed.
  Error: Feature write errors: Feature creation error (OGR error:
  Attempt to write non-linestring (POINT) geometry to ARC type
  shapefile.) Only 236 of 237 features written.

I think there is a geometry problem but I do not know how to isolate the problem or how to troubleshoot the error. Does anybody have any suggestions on how to fix this?

This is the same as the export to shapefile tool in ArcGIS


texts: Why vector tools cannot appear in QGIS data processing


texts: How to find vertices for a geometry in OpenLayers2

texts:
I have a LineString, from which I need to find the vertices logically and save those vertices. Can anyone help me out?

My code is:

OpenLayers.Geometry.LineString = OpenLayers.Class(OpenLayers.Geometry.Curve, { getVertices: function(nodes) {
    var vertices;
    if(nodes === true) {
        vertices = [
            this.components[0],
            this.components[this.components.length-1]
        ];
    } else if (nodes === false) {
        vertices = this.components.slice(1, this.components.length-1);
    } else {
        vertices = this.components.slice();
    }
    return vertices;
}
                }       );


texts: Passing variables between QGIS plugin and a combobox in layers custom form?

texts:
I am developing a plugin that asks the user to enter in "lineEdit" boxes, names of population types, depending on the number of types that they have previously defined.
Once these names are validated on the plugin, I want to be able to display them in a combobox within the form of the layer that was previously created by the plugin.
Basically, I want my plugin to be able to store its values and display them in a list within the form.
I have the code of my plugin in a file Chargementdialog.py and the initialization code of my form in shapefile_form_code.py.
Here is a part of the code I use in my plugin, in case the user chooses two types of Population:
class DeuxTypesPop(QtGui.QDialog):
  def __init__(self):
    QtGui.QDialog.__init__(self)
    self.ui = Ui_Dialog_Deux_Types_Population()
    self.ui.setupUi(self)

    self.ui.NamePop2_2.editingFinished.connect(self.store)

  def store(self):
    settings=QSettings()
    settings.setValue("/Nom_du_Plugin/config/Nb_Types_Population/Two/LPCD_One", self.ui.LPCDPop1_2.value())
    settings.setValue("/Nom_du_Plugin/config/Nb_Types_Population/Two/Name_One", self.ui.NamePop1_2.text())
    settings.setValue("/Nom_du_Plugin/config/Nb_Types_Population/Two/LPCD_Two", self.ui.LPCDPop2_2.value())
    settings.setValue("/Nom_du_Plugin/config/Nb_Types_Population/Two/Name_Two", self.ui.NamePop2_2.text())

And here is part of the code I used for the custom form
def shapefile_form_code(dialog,layerdid,featureid):
  global myDialog
  myDialog=dialog

  settings=QSettings()

  global PopType11
  PopType11=settings.value("/Nom_du_Plugin/config/Nb_Types_Population/One/Name", type=str)
  global PopType12
  PopType12=settings.value("/Nom_du_Plugin/config/Nb_Types_Population/Two/Name_One", type=str)
  global PopType22
  PopType22=settings.value("/Nom_du_Plugin/config/Nb_Types_Population/Two/Name_Two", type=str)

  global TypedePop
  TypedePop=dialog.findChild(QComboBox,"TypedePop")
  TypedePop.clear()
  NbTypePop=settings.value("/Nom_du_Plugin/config/Nb_Type_Population",type=int)

  settings.sync()

  global MyList

  if NbTypePop<7 and NbTypePop>0:
    if NbTypePop == 1:
      MyList=PopType11
      TypedePop.addItem(MyList)

    elif NbTypePop == 2:
      MyList=[]
      MyList.append(PopType12)
      MyList.append(PopType22)
      TypedePop.addItems(MyList)

When I put values in PopType12, PopType22, etc, instead of the settings.value... , the combobox gets updated, whereas by using this code, nothing happens, the combobox in the form stays empty.
I am new at Python and QGIS.


texts: Delete empty feature class in modelbuilder - "Calculate Value" Expression problems

texts:
I am trying to delete empty feature classes created in the process of intersecting layers with a set-up in modelbuilder.
By using Calculate Value tool i have written a code that should do the trick.

It is suppose to get count and if it is equal to zero then employ delete_management, which should delete the feature class.
Else, do nothing and add a warning.

Some way along, it is not doing what I want. The model is creating my feature classes but even though some of them are empty they are not deleted. The report states that:

" Start Time: Wed Mar 18 11:15:09 2015
WARNING 000117: Warning empty output generated.
Succeeded at Wed Mar 18 11:15:20 2015 (Elapsed Time: 11,44 seconds)
Executing (Calculate Value): CalculateValue DeleteIfZeroCount(r"M:\HAL\TeamNIS\Tools\SOE_check\Models\Version1\Results\TestIntersectOnWhole.gdb\NaturalFeaturesA_PortsAndServicesA_Intersect") "def DeleteIfZeroCount(inFC) :\n    import arcpy\n    if int(arcpy.GetCount_management(inFC).getOutput(0)) == 0:\n        arcpy.Delete_management(inFC)\n        return True\n    else:\n        arcpy.AddWarning('Error in '+inFC)\n        return False" Boolean
Start Time: Wed Mar 18 11:15:21 2015
The process did not execute because the precondition is false.
Succeeded at Wed Mar 18 11:15:21 2015 (Elapsed Time: 0,00 seconds) "

So, even though the feature class is empty, my Calculate Value tool doesn't delete the feature class.
I am a newbie at coding including arcpy, so I might have made an obvious mistake.
The question is why it won't delete the empty file? What is wrong with my code?
The model and code is shown in the attached image, please state if you need more info.


texts: Modifying polygon attribute interactively on web map for public interaction

texts:
I am looking for examples or ideas on how to create a web map using esri javascript api for redistricting purposes. On the web map, it would show the school boundaries in smaller sections (similar to census tracts) and would allow the user to interactively change the areas to create a new/modified boundary by clicking on the polygon and selecting another area using a dropdown list. It would then ask the user to print and/or email the new map to us for review. This map is the only thing I was able to find that is somewhat similar to what I am looking for: http://redistricting.greatergreaterwashington.org/

I thought about using the renderer unique value to separate my boundaries (ArcGIS API for JavaScript Sandbox) but I have not been able to find a way for a user to interactively be able to change the attribute to a different category.

Is this a possible task and does anyone have any examples I can use?


texts: Attributing length of line to first polygon but not overlapping polygons using PostGIS?

texts:
I have a polygon table that sometimes overlaps itself and has an integer column called rank

Then I have a line table

Basically I want to know what's the length of lines that are contained within each polygon. Problem is that because of the overlap, segments are being counted double for multiple polygons, making the lines look a lot lengthier than they actually are

Query used was

update polygontable a set length_of_lines =
(select sum(st_length(st_intersection(a.geom,b.geom)))
from linetable b where   st_intersects(b.geom, a.geom))

To avoid having lines look a lot lengthier, how to attribute the length of the lines according to the rank of the polygons?

Example: if a portion of the line is being accounted for polygon 1, it shouldn't be counted for polygon 2. In other words, if a line has 9 meters, 3 meters purely inside polygon 1, 3 meters in the area overlapping polygon 1 and 2, and 3 meters purely inside polygon 2 - How to update the table so that it attributes 6 meters to polygon1 and only 3 meters to polygon 2 (because the overlapping bit should be only accounted for polygon 1)


texts:
I'm trying to get satellite data from Hurricane Irma to asses damage. I tried the Landsat 8 sources but they don't produce tiles that can zoom close. I found better images from NOAA which are downloadable here:

https://storms.ngs.noaa.gov/storms/irma/index.html

My issue though is that there are many tif files that I need to get into Mapbox to create a layer. Is there a way to combine these and then upload them to Mapbox?

Here is the gdal info for one of the tifs.:

Driver: GTiff/GeoTIFF
Files: 20170911aC0813730w243600n.tif
Size is 18702, 18702
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-81.625100000000003,24.600100000000001)
Pixel Size = (0.000001347449471,-0.000001347449471)
Metadata:
  AREA_OR_POINT=Point
  TIFFTAG_DOCUMENTNAME=*
  TIFFTAG_MAXSAMPLEVALUE=0
  TIFFTAG_MINSAMPLEVALUE=1
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=0
  TIFFTAG_YRESOLUTION=0
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -81.6251000,  24.6001000) ( 81d37'30.36"W, 24d36' 0.36"N)
Lower Left  ( -81.6251000,  24.5749000) ( 81d37'30.36"W, 24d34'29.64"N)
Upper Right ( -81.5999000,  24.6001000) ( 81d35'59.64"W, 24d36' 0.36"N)
Lower Right ( -81.5999000,  24.5749000) ( 81d35'59.64"W, 24d34'29.64"N)
Center      ( -81.6125000,  24.5875000) ( 81d36'45.00"W, 24d35'15.00"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha


texts: Compress versioned database fails with underlying DBMS error?

texts:
We are after migrating 2 geodatabase from SQL server 2008, one to SQL Server 2014, the other to SQL Server 2016. Everything appeared to be normal at first but we just noticed that the weekly compress is failing on both databases. The compress process, run from ArcCatalog, runs for a few minutes and then fails with an underlying DBMS error. When I check the compress log I can see a record for that compress with a compress_start time but no finish time. The compress status is 'TRIMMING' and the before and after state count has only gone down by a small number.

Has anybody come across this problem before?


texts: QGIS: multiple mosaics bands to 1 Visual RGB

texts:
I have UAV raster converted with Pix4D mapper in transparents mosaics blue, green, red, nir, pir.
Which tool can I use to translate into a visual RGB image?


texts: Extending a linestring to the edge of an enclosing polygon in PostGIS

texts: I have linestrings that are each entirely within a polygon. I want to extend the first and last segments to touch the first rings of the polygon in line with them.


texts: Memory issue inserting a large dataset in GDB with ArcObjects

texts:
I am using the ArcObjects API for Java to load a very large dataset in a File Geodatabase in one feature class.

The whole dataset I have to load is composed of 200 datasets of 30k polygons, and I load the datasets one after the other. I read data from XML files (I am sure there is no memory leak in this step).

The loading process is in a BaseGeoprocessingTool.

Here is my code to write the data of one dataset in the feature class:

    IFeatureClass featureClass = ...
    IFeatureCursor insertCursor = featureClass.IFeatureClass_insert(false); // true uses too much memory
    IFeatureBuffer featureBuffer = featureClass.createFeatureBuffer();
    IGeometry geometry = null;
    List<Integer> fieldIdsToClear = new ArrayList<Integer>();

    for(FeatureBean featureBean : features) {
        if(log.canContinue()) {
            geometry = featureBean.geometry.getGeometry();
            geometry.setSpatialReferenceByRef(fromSpatialReference);
            if(toBeProjected) geometry.project(toSpatialReference);
            featureBuffer.setShapeByRef(geometry);

            // reset the fields in the featureBuffer
            for(Integer fieldId : fieldIdsToClear) {
                featureBuffer.setValue(fieldId, Variant.VT_NULL);
            }
            fieldIdsToClear.clear();

            featureBuffer.setValue(idFieldId, (int)featureBean.id);

            for(Entry<String, String> tag : featureBean.tagMap.entrySet()) {
                int fieldId = fieldIdMap.get(tag.getKey());
                featureBuffer.setValue(fieldId, tag.getValue());
                fieldIdsToClear.add(fieldId);
            }

            insertCursor.insertFeature(featureBuffer);
        }
    }

    insertCursor.flush();

    Cleaner.release(geometry);
    Cleaner.release(featureBuffer);
    Cleaner.release(insertCursor);
    Cleaner.release(featureClass);

After reading these advices, I use an insert cursor instead of IFeature.store() so that process is ~5x faster. So, I do not recreate the FeatureBuffer each time I insert a new feature. Instead of, I reuse it and clear the necessary fields.

This process works well and is very fast. Nevertheless, there is a memory leak. Watching ArcMap.exe process memory, we can see that it uses about 100 Mb to load each dataset, then the memory is not released, even if objects are released in my code. When memory is about 1300 Mb (only after loaded the 5th dataset), the process crashes on the exception AutomationException: One or more arguments are invalid

I have tried different Java max heap sizes with JavaConfigTool.exe between 512 and 1900 Mb, but it doesn't change anything.

What solution could I try to be sure that ArcGIS really release useless objects? As I have a lot of datasets to load, I really need to release to get back the 100 Mb each time a dataset has been loaded.

Any help is welcome to optimize my process.

EDIT: I use ArcGIS 10.2 with Background Geoprocessing (64bits) and have 8 Go of RAM.


texts: gdal_translate grib2 to geotiff raster bad values

texts:
I have used gdal_translate in the past with grib2 files without a problem, but when I run my command:

gdal_translate -of GTIFF -b 1 input_grib2 output_raster

My values for the output raster are -3 for all pixels (I believe this is the nodata value from the grib2 file).

When I bring the grib2 directly into ArcGIS I can see the actual values which range from -3 to 120.1

When I use gdalinfo -stats the values of band 1 (the only band) are all -3

I am trying not to use any arcgis (arcpy) commands in this script, and I am searching for a way to fix this just using gdal.

The grib2 file is coming from here:

http://mtarchive.geol.iastate.edu/2016/07/07/mrms/ncep/PrecipRate/PrecipRate_00.00_20160707-000200.grib2.gz

After clipping the grib2 file using wgrib2 -small_grib it works -- would like to find a way around this, but I needed a smaller area of interest anyways so it works out I guess.


texts: Displaying huge geotiffs (or vrts) with QGIS?

texts:
I am manipulating and processing global rasters at 30m resolution. The total raster sizes are usually [1,440,000 560,000]. I have access to a supercomputer, so I have written code that allows me to break up global rasters into manageable chunks, perform some calculations in parallel, and write them to disk quite rapidly.

I have hit a wall, though, when it comes to displaying results. I usually build a virtual raster of tiles covering the globe and pull it into QGIS. But it's incredibly slow (minutes to load, if it does). And if I try to pan or zoom, it's another many minutes. My first approach to solving this problem was to build overviews using gdaladdo. However, these take forever to build (as in days), which is not conducive to developing algorithms. Here is a list of the things I have tried and why/how they failed.

build overviews on the vrt. As mentioned above, this takes 2+ days to complete for 8 levels. That is unacceptable for my purposes.
build overviews on the individual tiles, then somehow merge into a vrt that contains the overviews. I am able to build overviews on the tiles rather quickly (supercomputer), but I have not been able to remerge them. I tried:

2a. gdal_merge on the tiles with overviews, but the overviews were not retained (or at least not recognized by QGIS) in the output tiff.

2b. gdalbuildvrt on the tiles with overviews, but as above, the overviews were not retained. [This is not correct, see edit.]

2c. I also tried a hybrid of building overviews for the tiles for levels 1-6 and building levels 7-8 directly on the vrt (basically option 2b) but it is still taking forever for just these two levels. I did some testing and see that the tile overviews are actually used to build to vrt overviews, but it's still on the order of a day to complete the overviews on the vrt.

So I am hoping someone here has some suggestions about where I should go next. Here are some options I'm considering:

Manually create the global pyramids myself. I am wary of recombining them into an .ovr file as I assume that will be tricky.
Use a mapserver (Geoserver). I know very little about this and am worried that it won't overcome the time hurdles while adding complexity to my process.
Split the domain by continents or some other region. I really want to avoid this option.

You might ask "why do you need to view the entire globe at 30m resolution?" One example: I take a mask of water pixels (globally) and skeletonize it to find rivers and perform measurements. My skeletonization algorithm requires a bit of tuning (for branch pruning, removing loops, general cleaning, etc.), and the output is necessarily at 30m. As rivers and landscapes are diverse across the globe, I need to be able to pan around to see the effects of any changes I've implemented.

I have also looked through QGIS to make sure there aren't any settings I could play with to render huge rasters faster, but I didn't see anything. Short of buying SSD drives, I think it's cranking as fast as possible. (My HDDs have I/O of ~250MB/s).

I discovered that building overviews on individual tiles, then building a vrt does apparently maintain the overviews--QGIS's "Pyramid" section in the metadata for the file is empty, but in the "Dimensions" section there is an entry for each level of overview (e.g. X 720000, Y 140; X 360000, Y 70, etc.). So I was wrong about 2b.

I also find that if I just pull all the tiles into QGIS, it renders in under a minute, while if I pull in the vrt that references the tiles, it takes >5 minutes (don't know how long exactly since I killed the process).

I did some testing on a computer with a SSD, and I found that I could load, display, and render the global vrts (without any overviews) sucessfully and at an acceptable rate. I have ordered a 1TB PCIe SSD in hopes that it will allow me to do the same on my computer. Will update with results.


texts: Map corners like asterisk symbol in QGIS?

texts:
I am trying to make a map like this image in QGIS.

I want to add min/max longitude and latitude and star symbol at the corners.


texts: Sourcing data from FME Server for ArcGIS Server?

texts:
I'm looking at potential architectures for a new project and the biggest unknown right now is how well ArcGIS Server can communicate with non-ESRI data sources for map services.

There are a number of data sources including Oracle Spatial and SDF files. There is no SDE and periodically converting source data into shapefiles or similar is not an option as the data must be live. It has been suggested that FME (server?) may be able to mediate between the two, but this is something I have never tried.

Does anyone have experience with this or a similar setup?

I'd be interested to know how stable a configuration this is, whether editing the source data would be possible (assuming a capable API was used on the client), or any other alternatives.


texts: Converting feature to raster, missing values

texts:
I have converted many features from the same shapefile to rasters, but am now getting an issue where values are missing. Converting the percentage values in column a212 (data type long) produces a raster with only 0, even though this feature has several rows with non-0 values. There are only 18 rows out of ~2.4 million with values higher than 0, but that shouldn't matter?

Having a similar issue with next column a213, where there are even more values >0, but is producing a raster with only 0 and 36?

See screenshot


texts: Leica Cyclone 5.1 TLS Data referencing using Trimble DGPS

texts:
I am using a Leica HDS3000 Terrestrial Laser Scanner (TLS) and Cyclone version 5.1 to acquire and process point clouds. Within the manual for Cyclone 5.1, Leica mentions that a user can enter the height of the TLS unit, the X, Y, and Z of a given point below the TLS unit.

Using a Tremble GeoXH, I can post process my location accuracy to 10 cm which is fine for my application. Given that I can get the X, Y, and Z from the Tremble for each station location; how can I project the TLS point cloud into a coordinate system such as WGS 84? It's not clear to me how this is possible within the documentation?


texts: Allowing different users to view and edit their workspaces with GeoServer web management interface

texts:
When different user loggeed into the GeoServer management interface, could they only view and edit different workspaces? Itried to do this by setting 'data security' and 'Users, Groups, and Roles', but it didn't work.

When I logged into the manager interface with the 'test' user, I got only the Layer Preview option in the Data section.

I have tried this answer ( How to allow users to publish data in their workspace? ), but still have this problem.

How should I solve it?


texts:
I need support on displaying the form (designed by using QT designer) on top of the QGIS window.

My concern is, after loading my plugin, my designed form is opening on top of the QGIS application window. If i click any other place on QGIS window, the opened UI window disappears and moves behind QGIS. I want to access UI window and QGIS canvas same time.

How can I design a UI like this?


texts: Calculating Driving Times from Cities to each Raster Cell in State using ArcGIS Network Analyst?

texts:
I am doing a project for a depart of the Utah state in ArcMap 10.4.1.

The project is create a model that will value each raster cell for the state of Utah with driving times from each of the Top 2 Populated cities in each county. I have the county layer, but am having a really hard time trying to get the Network Analyst work.

Any ideas where to start?

I have tried to do the service area in the Network Analyst, and load the Facilities as the cities layer, but it keeps saying the "cities" aren't found in the Network Analyst layer. Additionally, I am not sure how to have it run driving times from each city to each raster cell. Most of the Network Analyst work has one or a few desired destinations. After I have all the cells valued for each driving time from each city, I will just take the max value.


texts: Why does my GeoJSON have no coordinates when saving from Shapefile in QGIS?

texts:
I've got a shapefile (.shp, .shx, .dbf, and .prj) that I'm trying to save as a GeoJSON using QGIS 2.14. When I load the .shp into QGIS, everything seems to work fine, the Shapefile renders properly.

However, when I try to do a Save As and select GeoJSON, I get null values for all my GeoJSON coordinate values.

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { "Section": null, "Block": "Last Name, First Name" }, "geometry": { "type": "MultiPolygon", "coordinates": null } },
{ "type": "Feature", "properties": { "Section": null, "Block": "Last Name, First Name" }, "geometry": { "type": "MultiPolygon", "coordinates": null } },
{ "type": "Feature", "properties": { "Section": null, "Block": "Last Name, First Name" }, "geometry": { "type": "MultiPolygon", "coordinates": null } },
{ "type": "Feature", "properties": { "Section": null, "Block": "Last Name, First Name" }, "geometry": { "type": "MultiPolygon", "coordinates": null } },
{ "type": "Feature", "properties": { "Section": null, "Block": "Last Name, First Name" }, "geometry": { "type": "MultiPolygon", "coordinates": null } },
etc, etc..............

Is there something special I need to do to make sure my coordinates show up when I'm saving my data?


texts: Overview effect on postgis polygon table

texts:
Is There a Way, with postgis and qgis, to reproduce the overview effect of tiff Raster on a polygon table in order to increase speed?
What if we could create simplified versions of this table,gathered under the same view.
Until now, i only play with different scale visibility on different simplified version of the table.


texts: What is the logic behind the order of features when getting features by PyQGIS?

texts:
When I run the script below, I get features as ordered by id as expected.

layer = iface.activeLayer()

# getFeatures()
ids = [f.id() for f in layer.getFeatures()]
print(ids)

# OUTPUT: [0L, 1L, ..... 499L, 500L]

If I select all features and run the next script, I get the same output.

# selectedFeatures()
ids = [f.id() for f in layer.selectedFeatures()]
print(ids)

# OUTPUT: [0L, 1L, ..... 499L, 500L]

But if I select manually some features -for example features with 12, 13, 17 id in my shapefile- and run the previous script, I get unordered list of id. I always get the same output, no matter what the selection order is.

ids = [f.id() for f in layer.selectedFeatures()]
print(ids)

# OUTPUT: [17L, 12L, 13L]

## OR selectedFeaturesIterator()
ids = [f.id() for f in layer.selectedFeaturesIterator()]
print(ids)

# OUTPUT: [17L, 12L, 13L]

Why selectedFeatures() and selectedFeaturesIterator() return unordered features by id, while getFeatures() returns ordered features? What is the logic behind?

EDIT: The point I wonder is not whether features are ordered in database/file/memory or not. What I mean by "why" is that why a method returns ordered features but another one returns unordered or same method returns different results. Pay attention to selectedFeatures(). As an example, if all features are selected, it returns [1, 2, 3, 4]. But it returns [2, 1, 3] if some features are selected.

NOTE:

getFeatures() and selectedFeaturesIterator() return an iterator (QgsFeatureIterator), selectedFeatures() returns a list (QgsFeatureList).
QGIS 3 has getSelectedFeatures() instead of selectedFeaturesIterator() in QGIS 2.
I get the same result in both QGIS v2.18 and v3.2.


texts: QGIS 3 - load point layer from Geopackage doesn't render points

texts:
I am trying to load a point layer in QGIS (3.0.2) using DB Manager. I have a non-spatial table, which I imported to the Geopackage from CSV (and needed to implement the .csvt approach to import the column types correctly - https://anitagraser.com/2011/03/07/how-to-specify-data-types-of-csv-columns-for-use-in-qgis/ )

Using SQL, I have worked out how to create a geometry column from two in my table (lat, long) using makepoint(), which means I should be able to create a point layer with a subset of my table columns and load in QGIS:

Select Mykey, makepoint( "P1_location_Longitude" , "P1_location_Latitude") as geometry, "Photo1" from MyTable

However, when I load the layer it does not render any points (even though the attribute table has all values from the parent table). When checking the layer properties the symbology tab is blank (i.e. no symbol rendering is possible) and the Error console shows the message no renderer for this layer.

Same behaviour in QGIS 2.18.3.

How do I specify a symbology or allow the points to be rendered for this layer?

If I create an Sqlite database, and complete the same operations (i.e. import the CSV, use the same query and the same settings in DB Manager) the layer is rendered as expected, so I wonder if it is an issue with Geopackage format in QGIS?

I haven't found any reports of this being a problem / bug so I'm assuming that I'm missing something.


texts: Why doesn't QGIS show the results from a WFS BBOX from a Python script

texts:
I'm trying to figure out how to do get requests in QGIS using a BBOX. I can get a BBOX to clip a layer when I paste my URI into a browser but the EXACT same URI returns an empty layer in QGIS.

A simple easy to try specific example: The following URI returns four features when pasted into my browser.

http://ogi.state.ok.us/geoserver/wfs?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=ogi:quad100&BBOX=-100,37,-99,38,EPSG:4326&SRSNAME=EPSG:4326

However the following lines create an empty layer (either pasted in to the console or run as a script)

sourceUri = "http://ogi.state.ok.us/geoserver/wfs?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=ogi:quad100&BBOX=-100,37,-99,38,EPSG:4326&SRSNAME=EPSG:4326"
vlayer = QgsVectorLayer(sourceUri, "QUADS", "WFS")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

Interesting to note that without the crsurl I need to swap the lat/lon order to correctly clip data in the browser. with or without the crsurl and which ever lat/lon order I use in QGIS I still get an empty layer.

As far as I can tell my URIs are correct (they certainly generate a correct result in a browser), so why doesn't QGIS get any data?

It would be a shame to have to write my own HTML/XML post/parse script in Python just because I'm missing some small aspect of QGIS.


texts: How to 'fill' a region in a raster layer specified by a shape with a value

texts:
I have a shapefile for a country, made up of a shape per province.  I want to create a raster layer in the shape of the country and for each pixel in the area covered by province 1 set the value to value1, for the area covered by province 2, set the pixels to value2, etc.  I want the pixels to represent 1km/1km resolution.

I've tried various ways to go about this, but the closest I've come is some combination of creating a constant raster layer and cropping it to the whole country boundary, and I've had some trouble getting anywhere with that and suspect it's not the best way to go about it.

I'm using QGIS, but can switch to anything open-source that helps (GDAL, R, etc).


texts: GeoServer 2.18. Setting Min & Max Scales

texts: I'm using CSS and wanted to know how to set the min and max scale for a given layer.  I see many solutions referencing SLD, but want to do this in CSS.  I really need help with the syntax.


texts: Calculate Grid Convergence Angle within ModelBuilder Model

texts:
I am trying to use the 'Calculate Grid Convergence Angle" tool (from Cartography toolbox -> Data Driven Pages) into ModelBuilder. The tool is pretty simple -

(1) import,

(2) remove no data

(3) project to LCC and then

(4) recalculate grid convergence angle based on a field.

I keep getting an error

  error 000229 [The data is not accessible. It is likely either invalid
  or not available]

Why is it unable to use the output from step (3)'Project to LCC output' (which is populated with the necessary fields in the attribute table) as an input for step (4)?

Edit: image 2 of CalcGCA tool added.

edit: more images added:


texts: Is it possible to create a data frame graphic using arcpy?

texts:
I'm interacting with an MXD using arcpy, and wish to create a graphic within a data frame - as with the Draw toolbar when working manually. The GraphicElement page says...

  The GraphicElement object provides access to properties that enables
  its repositioning on the page layout as well as methods that allow for
  duplicating and deleting existing graphic elements.

...which doesn't mention the ability to create a new graphic element. Is there another area I should be looking at in arcpy, or is this not currently supported?


texts: How do i select specific time and download quality MOD021KM data

texts: I am trying to  download MOD021KM for net primary production estimation, and am so stuck on how to download the data and selecting the best time of a day using LAADS website. Morever which values for tile or longitude/latitude should i enter to obtain accurate geographical location for Zimbabwe?


texts: Find closest points using S2, given point and radius

texts:
I'm trying to use the go s2 package to find closest points, given a point and a radius (meters).  I'm modeling my alg. after dynamodb-geo.js as I'm storing my points in DyanmoDB.

I can't seem to get a non-empty list of covering cells (coveringCells := rc.Covering(r) below).

What am I doing wrong in my algorithm? Is there an easier way to achieve this?

I think I could use RectFromCenterSize() to compute the rect given a center lat/lng and a radius - but I'm not quite sure how to compute the size param for this method.

// Query a circular area constructed by a center point and its radius.
// @see https://github.com/rh389/dynamodb-geo.js/blob/6c388b9070014a096885e00fff6c3fc933d9853f/src/GeoDataManager.ts#L229
func queryRadius(lat float64, lng float64, radiusMeters float64) (error) {
    earthRadiusMeters := 6367000.0

    // Step1: Get the bounding region (rectangle) from the center and the radius
    // @see https://github.com/rh389/dynamodb-geo.js/blob/6c388b9070014a096885e00fff6c3fc933d9853f/src/s2/S2Util.ts#L23
    centerLatLng := s2.LatLngFromDegrees(lat, lng)

    latReferenceUnit := 1.0
    if lat > 0.0 {
        latReferenceUnit = -1.0
    }
    latReferenceLatLng := s2.LatLngFromDegrees(lat+latReferenceUnit, lng)

    lngReferenceUnit := 1.0
    if lng > 0.0 {
        lngReferenceUnit = -1.0
    }
    lngReferenceLatLng := s2.LatLngFromDegrees(lat, lng+lngReferenceUnit)

    latForRadius := radiusMeters / centerLatLng.Distance(latReferenceLatLng).Radians() * earthRadiusMeters
    lngForRadius := radiusMeters / centerLatLng.Distance(lngReferenceLatLng).Radians() * earthRadiusMeters

    minLatLng := s2.LatLngFromDegrees(lat-latForRadius, lng-lngForRadius)
    maxLatLng := s2.LatLngFromDegrees(lat+latForRadius, lng+lngForRadius)

    boundingRect := s2.RectFromLatLng(minLatLng)
    boundingRect = boundingRect.AddPoint(maxLatLng)

    // Step2: Compute the CellIDs for the region we want to cover.
    // defaults per https://github.com/vekexasia/nodes2-ts/blob/1952d8c1f6cb4a862731ace2d5f74d472ec22e55/src/S2RegionCoverer.ts#L101
    rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 8, LevelMod: 1}
    r := s2.Region(boundingRect.CapBound())
    coveringCells := rc.Covering(r)

    for _, c := range coveringCells {
        log.WithFields(log.Fields{
            "Covering Cell": c,
        }).Info("=>")
    }

    return nil
}


texts: Concatenate fields (double and string) to create address using arcpy

texts:
fields:

STNUM = double
STRDIR = string
STRNAME = string
STRTYPE = string

Tuple and s.format

I have tried:

arcpy.CalculateField_management (newparcelsUTM, "SiteAddress","{0} {1} {2} {3} ".format('!STRNUM!','!STRDIR!','!STRNAME!','!STRTYPE!'),"PYTHON_9.3")
## set variable- str(double) ##

I have tried setting a variable to change the double field to a string:

AddNum = str('!STNUM!')

then use it in the Calculate field line:

arcpy.CalculateField_management (newparcelsUTM, "SiteAddress", AddNum+" "+'!STRDIR!'+" "+'!STRNAME!'+" "+'!STRTYPE!'),"PYTHON_9.3")

All my attempts end with:

  ExecuteError: ERROR 000539: SyntaxError: invalid syntax (, line 1)
  Failed to execute (CalculateField).

I have tried different iterations on this- single quotes v double quotes... it works if I use just the STNUM field. It works if I use any combination of the string fields. It just does not work combining the double field with the string fields! I have read many threads regarding concatenate, Field calculator, combine double with string, on stack exchange and the Esri sites... I have been trying!


texts: Performing raster analysis using Arcpy - in_memory does not preserve raster shape

texts:
I am creating a script tool using Arcpy and it transforms elevation data into a final clip polygon as part of a lengthy site selection tool. When I perform the raster analysis and save my intermediate data to disc, it works fine; however, when saving the intermediate layers to in_memory, it does not preserve the shape of the output rasters. It forces them into a rectangle.

Slope_1 = "D:\egm722\Data\Scratch.gdb\\Slope_1"
arcpy.gp.Slope_sa(Elevation_Data, Slope_1, "DEGREE", "1", "PLANAR", "METER")

Slope_1 = "in_memory\\Slope_1"
arcpy.gp.Slope_sa(Elevation_Data, Slope_1, "DEGREE", "1", "PLANAR", "METER")

How can I get the shape preserved using in_memory? This is only one process of the tool which contains many processes and I don't want my end user to have to save lots of intermediate data to disc.


texts: Point to Raster Arcmap - incorrect pixelvalues?

texts:
I am having a problem using the arcmap point to raster tool. I have x/y coordinates for a set of data, and 1464 (2 months worth of hourly data) separate files to convert from an excel with X/Y/value to a matrix for use in another program. I already have a small python-script in order to automate the point to raster and then a raster to ascii, but in my first round of testing, the values i get out of this are really weird.

I have 6 separate excel files with x/y/value1/value2/value3... etc., which i have loaded into arcmap using the Add Data -> Add xy Data tool. However, when i convert this point data based on a value-field, it gives me a range of 0:401, whereas the actual maximum value in excel is 5.05. How can these values be so different?


texts: Running a new iteration for each column of a dataset in ArcGIS

texts:
I'm pretty new to Modelbuilding, having just been turned onto it about a week ago on here. So my problem is, I currently have a dataset of habitats that change through time (from Immature to Mature, or remaining Mature, etc). I would like to create a shapefile for each habitat type for each year. Currently I am able to do this by hand, running a Feature Selection iteration on each field and running the "feature class to feature class" tool to export it as a shapefile. What I was wondering was, is there any way to automate this process across a number of different columns? So, run the iteration and shapefile export for the 1995 field, 1996 field, 1997 field, etc, etc. I'm trying my best right now to sort it out in Python, but today was day 1 of learning and was hoping that there was another method out there that someone may be able to point me in the right direction of.

Hopefully this makes sense, if you need any more info, please let me know.

Thanks very much in advance.

EDIT - Example data - The cut off for going from young (or Disturbance) to a mature forest is 40 years. So the field "Year" is the year the forest began growing (or year since it was removed) and under each year is the forest class for that year.

FID    Year    1995           1996          1997         1998
0      1992  Disturbance  Disturbance  Disturbance  Disturbance
1      1854  Black Spruce Black Spruce Black Spruce Black Spruce
2      1960  Disturbance  Disturbance  Disturbance  Disturbance
3      1875  Deciduous    Deciduous    Deciduous    Deciduous
4      1956  Disturbance  Disturbance  Jack Pine    Jack Pine
5      1957  Disturbance  Disturbance  Disturbance  Lowland Conifer

And this is the model I am currently using, I would like this to be re-run across all of my year fields (columns)


texts: How to draw polygons in QGIS perfectly aligned

texts:
I would like to draw polygons for an indoor map test and would like to draw perfect rectangles (horiz. and  vertical edges).
How to set QGIS options on polygon edit /draw to do that?


texts:
I am trying to set up a wms server using point data in a postgis (v 1.5.1, pgsql v 8.4)  database and mapserver (v 1.5.1) (I know these are not the most current versions of anything, but my server is running debian stable, if neccesary, I'll upgrade them).

The data is a global data set, so the points are stored in lat/lon as postgis geography data type. I manage to display the data (using either arcmap 10 or Qgis 1.8 or directly in a web server) in lat,lon, but I would like to be able to serve the data in other projections as well. I have tried to set up LAEA north pole (EPSG 3408) as well. In the GetCapabilities, I get a part like this:

<Layer>
    <Name>xxxx</Name>
    <Title>yyyy</Title>
    <Abstract>zzzz</Abstract>
    <CRS>EPSG:4326</CRS>
    <CRS>EPSG:3408</CRS>
    <EX_GeographicBoundingBox>
       <westBoundLongitude>-50</westBoundLongitude>
       <eastBoundLongitude>90</eastBoundLongitude>
       <southBoundLatitude>-50</southBoundLatitude>
       <northBoundLatitude>90</northBoundLatitude>
    </EX_GeographicBoundingBox>
    <BoundingBox CRS="EPSG:4326" minx="-50" miny="-50" maxx="90" maxy="90" />
</Layer>

I.e my two projections are listed, but I only get the bounding box for the original projection. SHould I also get the other one?

As far as I can see from the instructions I can find, I just need to list the projections I want to use in wms_srs in the METADATA section - in the layer, in the MAP or both places (I've put it both places, but I've also tried it in just the layer or MAP, with the same result)


texts: Switch to another grapic card on my Laptop for ArcGIS

texts: I use ArcMap on a laptop, a Lenovo Yoga P40. This laptop has two graphic cards: the standard Intel graphic card and a NVidia graphic card designed to be used on graphic softwares (e.g. AutoCAD Inventor). Currently, I think my laptop runs ArcMap with the Intel graphic card, but sometimes I notice it has issues to load and handle ArcMap, therefore I was wondering: is there a way to set my laptop to run ArcMap with the other graphic card?


texts: Arcpy Dissolve creating grid that is not produced when run manually via ArcMap

texts:
Software: ArcMap 10.5.1, ArcPy

I have created 2 python scripts that seeks to create a floodplain polygon from a DEM. py1 is DEM to Polygon while py2 cleans polygon to create final output. In py2 there a number of intermediatory steps to produce and clean the polygon to get the final floodplain product (clip, dissolve, feature removal, eliminate, etc.). Dissolving the Floodplain polygon by a water catchment to get single part features is one of these steps.

Before creating the script, I always run the steps manually in ArcMap to ensure processes/outputs are correct and steps are working. Following this, I generate code snippets and populate a working script. Everything is running and producing a final output.

However, I have noticed that when I review the output for the dissolve I am getting a kind of "gridded" output that I do not get when ran manually in ArcMap (either by using python console or the tool itself). This obviously has an affect on subsequent steps as this grid is visible for each step following the dissolve.

Image below shows dissolve outputs of ArcMap (1 feature selected) vs ArcPy (2 features selected). Notice the obvious right angles from the "grid"

I have checked environments for both ArcMap and ArcPy process and they are pretty much the same (apart from workspaces. I have added in additional settings to script to fully match those in ArcMap [outputCoordinateSystem and geographicTransformations]):

Code snippet:

#2 Dissolve by gridcode to remove multipart features
input_p2 = output_p1
output_p2 = os.path.join(output_root + "_step2_diss")
diss_field = "gridcode"
print ("Dissolving '%s'..." % (os.path.basename(input_p2)))
arcpy.Dissolve_management(input_p2, output_p2, diss_field, "", "SINGLE_PART", "DISSOLVE_LINES")
print ("Created: %s" % (os.path.basename(output_p2)))

link to dissolve tool - http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/dissolve.htm

This is not isolated to my PC as colleagues and the client have tested on their PCs both in and outside of our Company. Has anyone encountered this issues before?


texts: Ways of editing Postgis tables from remote desktops (WFS?)

texts:
I've got like 500 postgis tables and i would like to serve then as WMS and especially as WFS to be able to edit that info. To do so i read that the easiest way is Mapserver + TinyOWS. Mapserver works with .map files so i need to generate the map files form postgis tables.

I am also testing geoserver. But i am getting a Provider errors: empty response whenever I try to edit the layer.

Which would be the easiest way to make a fully working WMS + WFS service?

Regards,


texts: Using QGIS to sort datapoints into districts

texts:
I have a dataset of about 15000 events (protests) across the US with information about which state, county and city each event occurred in.
The problem is I also need to know which congressional district each event occurred in. Neither county or city on their own can reliably tell me about congr. district, but the dataset does however also contain coordinates (lat/long) for each event.
If I've understood things correctly I need to find shapefiles (?) for US congr. district (which I have found), extract some information from them about a geographic definition (?) of each district, which I can then load into R (the data software I'm using) and then use R to do the assigning.
I also understand that I need to make sure the coordinates in my dataset use the same coordinate system (coordinate reference system?) as the shapefiles (or something like that).
Have I understood the basics of what I need to do right? Is it a very complicated thing to do? Could anyone help me out or point me in the right direction?


texts: Looping over array of type float with interval of float (NDVI) using Python?

texts:
I am stuck with my Python Script. The case is, I have an array of floating numbers from 0 to 1 (NDVI values). I want to loop over my array with an interval/step of every 0.025. Every interval/step I have to calculate something a+(b*element in the loop).

vNdvi = one dimensional array of NDVI values coming from a MODIS image

rNdvi = vNdvi[np.logical_and(vNdvi>=0.3, vNdvi<0.8)]

rNdvi.sort()


texts: To import bulk data in google fusion table

texts:
I want to import bulk data into fusion table . What is the best data format should I use to import data i.e JSON ,CSV .

If I want to import data in JSON format how can I do it since on google help only way was to upload CSV files .


texts: Create new regions in layer by filling in between polylines on existing layer

texts:
I'm not sure if this is possible but would save me hours of work.

I have a layer which consists of polylines. Most of those polylines bound areas; houses, roads, etc.

Currently I am required to colour each of those areas to define what they represent in the real world. My current technique is drawing a region in a new layer; as you can guess this takes a lot of time, especially for complicated areas.

Is there a feature in QGIS whereby I can create a new region on a layer by selecting a point within the bounded area and it creates the region shape?

EDIT

I have tried converting the whole layer of polylines to regions but the results were not as desired, there are some lines that are not bound to an area (partial fences, walls etc) so the results look odd.

It would be something similar to 'hatching' in an CAD program, to fill the area bound by lines. Only I want to do this to create the region.


texts: Projecting from Lambert azimuthal equal area (LAEA) projection to WGS84 without distortion

texts:
I need to convert a shapefile from Lambert Azimuthal Equal Area projection system (North Pole) to WGS84. When I use QGIS to do this, the lines that cross the anti-meridian cross over (horizontal black lines in the image).
To avoid this, I found one solution that worked, which is clipping the region of interest and making sure the region I clipped does not include the -180 and 180 longitudes (the reason for distortion apparently), but is there an out-of-the-box tool (say from some standard Desktop GIS or GDAL or so) to deal with this problem without doing manual clipping or having to program it myself?

Source projection (Lambert Azimuthal Equal Area projection):

Target projection (WGS84) that looks distorted:
.

Solution that worked for the region of interest:

The '-wrapdateline' option in ogr results in a similar but slightly different distorted output.


texts:
I would like to convert all values having the value 1 in a binary raster to a line vector layer. The following image shows the binary raster (1=white and 0=black):

The optimal output would look like the red lines in the following image (which I realized by modifying the GDAL script raster to vector line):

However, this is very slow and not applicable on large rasters. Is there a similar open-source implementation (callable in Python) that does the same with less processing time?
Would it decrease the processing time of the GDAL script raster to vector line when writing each single wkbLineString to a shapefile instead of adding them all to a wkbMultiLineString before writing them to a shapefile?
The only related posts I found so far are referring to polygonizing the raster (which I realized with GRASS GIS r.to.vect), which creates polygons and thus two lines, where I need only one:


texts: Creating simplified (binned) polygons from raster file?

texts:
I have esri ascii grid raster files like

ncols 4
nrows 3
xllcorner 0
yllcorner 0
cellsize 1
nodata_value -999
1 2 3 4
4 21 3 3
3 2 31 1

How can I convert this into vectorized (WKT) format when simplifying / reducing polygons via a binning of the values i.e. neighbors within a range of 0...10...20...30 should be put into the same polygon.

Naively this translates to

1.0;POLYGON ((1 3, 0 3, 0 2, 1 2, 1 3))
4.0;POLYGON ((1 1, 1 2, 0 2, 0 1, 1 1))
3.0;POLYGON ((0 1, 0 0, 1 0, 1 1, 0 1))
2.0;POLYGON ((1 2, 2 2, 2 3, 1 3, 1 2))
21.0;POLYGON ((1 1, 2 1, 2 2, 1 2, 1 1))
2.0;POLYGON ((1 1, 1 0, 2 0, 2 1, 1 1))
3.0;POLYGON ((3 1, 4 1, 4 2, 3 2, 3 3, 2 3, 2 1, 3 1))
31.0;POLYGON ((3 1, 2 1, 2 0, 3 0, 3 1))
4.0;POLYGON ((3 2, 4 2, 4 3, 3 3, 3 2))
1.0;POLYGON ((4 0, 4 1, 3 1, 3 0, 4 0))

I can aggregate/bin the numbers into a range from 0...100 with steps of 10 and group-by aggregate (ST_CONCAT/ST_UNION) the geometries, but this is not what I want to achieve. I.e. When creating the polygons I would rather only want to create them in aggregated fashion. So for the value range, it should be already in the binned range, and for the polygons they should be considered neighbors if a) next to each other b) binning interval is the same.


texts: Aggregating elements of array using QGIS

texts:
In QGIS 3, are aggregate functions over a set of array elemnts possible?
As simple as in the form of:
aggregate(array:=array(1,2,3), aggregate:='mean')
I'd use it for the task to get an average for an array of raster values, which I extract at the nodes of a linestring input feature layer. The raster values at the linestring nodes I can retrieve beautifully with the below expression, but obviuosly, there is no way to aggregate those values into on single aggregate for each feature.
with_variable('nodes',
    nodes_to_points( $geometry),
    array_foreach(  generate_series(1, num_points(@nodes)), raster_value(DGM_Slope', 1, point_n(@nodes, @element)))
)


texts: Import TomTom OV2 into Google Maps as My Place Stars

texts:
Is there anyway to import my TomTom POVs (saved in an OV2 file) into my google maps account as favorites? i.e. My Places > Starred.

I suspect this may involve converting the OV2 into something else and possibly using google earth or other software.


texts:
I want to join 7 images [that is in .tiff] and merge and than subset the area using shapefile. But i am not able to merge it. I have loaded the images as raster layer and than merge it Using Qgis and not getting desired result.
This is the sample image when loaded as a raster layer.

Adter I load all the iamges it looks like this.

So which software which method should i used?

I am getting this result and the color is blown gone out.


texts:
We have a table in MariaDB that contains a lot of point of interest. These POI have a list of numerical attribute such has occurrence, weight, etc. We are looking for a way to aggregate these POI by tile (bounding boxes) to finally display the result on a map (google maps).

Is it possible to do so in a request like : SELECT SUM(...) GROUP BY bounding_boxes?

*Note that we can use PostGIS if there's no way to accomplish these kind of request with MariaDB. (but I don't even know if it is possible)


texts: ArcGIS 10.2.1 working with time and reporting when feature geometry has changed

texts:
I'm working, with ArcGIS 10.2, ona whole lot of petroleum licencing data that chops and changes, in active areas, on a relatively frequent period. i collect the geometry data month by month (saving the copy as a feature class) and then copy the new geometry outs as new features to a feature class that represents the accumulation of geometries over time.

The product, described below, is designed to answer: When has this licence's geometry changed? This is important as companies change their licences for a reason either through business deals, partnerships or failed ventures etc...

I have created a product that will show 90-captures of a geometry due to the time period data has been gathered for. However, due to the nature of the licencing data a lot of licence areas remain unchanged for long periods of time thus making working with the data a bit cumbersome.

For now it is manageable for me to use definition queries to sample a single licence and with the time slider play the history stopping when the geometry changes, making a note of this for later for further analysis (i.e. presenting the time stamps for change as well as a map of the geometry change)

So, after setting the scene the question is: Is there a way, using Python, to print the date (an attribute field) of a licence for when its geometry changes?

or

Is it possible to select a feature in a time-enabled layer and detect when its geometry changes?


texts: Arc error executing function Skyline with surface

texts:
Trying to get a skyline for a set of 3D points based on a DEM (LiDAR based raster), but get the catch all error 999999.
The surface is optional and the function works when no surface is assigned. All other optional settings are set to default.
I've found the same problem posed elsewhere, but no solution yet (GeoNet and Shadow Analysis Tool Comments).
I've attempted in Arcmap 10.1 and 10.2.2. I've also tried copying the raster to into a .gbd, converting the raster to integer, and executing the function from python with no luck. Am I missing something obvious?


texts: Maintaining Borders when Merging shapefiles?

texts:
I am trying to merge two layers of a map which includes the state and county levels. I was able to edit the borders to look how I want, but when I merge them I am struggling to keep the style formatting.

This image is how I want them to look, with State level laid on top of Counties. I would like to export this as one shapefile.

When I merge the Shapefiles, the map ends up like this. When I edit the borders of the states it now changes all the borders to thick black.

Is there a way I can keep the formatting from how Image 1 looks, while merging the two shapefiles.  Below are the two shape files without being overlaid.


texts: I'm working on a dataset from an excavation collected by an ancient total station. I have a point layer with point code, xyz values and a description of the point. These points are corners of features, lines around features etc etc. They need to be converted to polygons/lines to make sense of all the different features. The problem is that this site is extremely complicated because it consists of a lot of (partial) inhumations on top of eachother. A 3d visualisation of this data would be perfect to analyse this dataset. A 3d visualisation of the points is easy using qgis2threejs. I would like to know how i add the z values from the points and add them to (to be created) polygons so i can create a 3d visualisation of the polygons?


texts: How to create polygon that goes over +/-180 (dateline) in ArcMap?

texts:
I want to create a polygon feature that will go through the dateline (180 long.) in ArcMap. I am using WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857) for my data frame and my polygonal feature class is defined with the same projected coordinate system. My feature classes are stored within enterprise geodatabase 10.2.1 (SQL Server Geometry spatial type). However, the same behaviour can be reproduced with the file geodatabase feature class.
From Esri help:

Data in a geographic coordinate system is treated, in a sense, as
two-dimensional, at least to the point that there are edges, just as
in a projected coordinate system. Generally, the edges are at
-180/+180 for the east–west extents, and -/+90 for the north–south extents. These extents are called the "horizon" of the coordinate
system.
In a geodatabase, a feature that crosses an edge will be clipped to
the horizon. If you are editing features in ArcMap, the data frame's
coordinate system will affect what happens to a feature that crosses
the +/-180° line. If the data frame is using a geographic coordinate
system, that is, the +180 and -180 lines are separated, a feature that
crosses outside the horizon will be clipped to the horizon. If the
data frame is using a projected coordinate system, with the +/-180°
line inside, a feature that crosses it will be split into multiple
parts.

I am able to create a polyline feature that crosses the 180° line. A multipart feature is created just as outlined in the Help. When I run the Multipart To Singlepart GP tool, I can verify that my line was originally multipart and is represented in the output feature class as two singlepart features.
However, this doesn't happen when I am trying to create a polygon feature crossing the 180° line. I am getting the error message:

Create
The Create Feature task could not be completed. The coordinates or measures are out of bounds.

Why don't polygons get splitted just lines did? Is there any workaround that would help me to achieve my goal?


texts: Geoserver SLD: link to symbol in database instead of file

texts:
I have a PostgreSQL table with fields:

'id' with identificators of objects types (e.g. shops, cinemas etc),
'graphic' with SVG codes of correcponding map symbols.

I want to create sld style with links to graphic in this table, instead any files.

EXAMPLE

Classic style:

<PointSybmolizer>
    ...
    <Graphic>
        <ExternalGraphic>
            <OnlineResource xlink:type="simple" xlink:href="sign.svg"/>
            <Format>image/svg+xml</Format>
        </ExternalGraphic>
    </Graphic>
</PointSymbolizer>

And style I want:

<PointSybmolizer>
    ...
    <Graphic>
        <!-- Link to database record with same sign.svg -->
    </Graphic>
<PointSymbolizer>

Are there methods to achieve this?


texts:
I am trying to clean existing items in the rubberband. The codes work fine on QGIS 3.2 but on 3.10 crashes.
rbs = [ i for i in canvas.scene().items()
                if issubclass(type(i), qgis._gui.QgsRubberBand) ]
if len(rbs)>0:
                for rb in rbs:
                    if rb in canvas.scene().items():
                        canvas.scene().removeItem(rb)

Can someone explain what is wrong with this code?
I have another function drawing on the lines on the canvas and this function should remove those.
To regenerate my issue you can run the following code.
 #to draw polyline
canvas=iface.mapCanvas()
r = QgsRubberBand(canvas, False)  # False = not a polygon
points = [QgsPoint(-100, 45), QgsPoint(10, 60), QgsPoint(120, 45)]
r.setToGeometry(QgsGeometry.fromPolyline(points), None)

#the next lines are the part I have problem with
rbs = [ i for i in iface.mapCanvas().scene().items()  if issubclass(type(i), qgis._gui.QgsRubberBand) ]
for rb in rbs:
        if rb in iface.mapCanvas().scene().items():
            iface.mapCanvas().scene().removeItem(rb) # this line causes QGIS to Crash

I have tested this code on several PC and QGIS. It works perfect on every PC upto QGIS 3.2 and it causes a crash for every version after.


texts: How to match GPS collected data to routes on Google Maps?

texts:
I am using Google Maps to implement a GPS Monitoring system of data collected from vehicles. My problem occurs when the points and lines do not match with the route and cut through other streets and structures such as shown in the example below.

How can I solve this problem?


texts: How to add a column to "link" lines to polygons (Qgis)?

texts:
First i'll show you this image

The green polygons are lots/land/property of different people of my city.
The blue lines are the river banks of a river that crosses my city.
The red lines are called here as "Estacas"...  in english they can be translated as "poles".

What are these poles. They are lines that  measure the distance of a construction (like a house) to the river bank.

What I have to do now is link, somehow, these poles to the lots.... each lot have an ID and each pole also have an ID. But there is lots of poles per lot. Like, the first lot, right-to-left, has seven poles crossing it. So, these seven poles, and only them, have to be linked to that one lot. (and so on to with the other lots and other rivers on my city)

Can I do that?
There's a way, with Qgis, to do it?


texts: QGIS not loads a layer from Geoserver WFS

texts:
This is my geoserver http://webservice.dftrans.df.gov.br:8080/geoserver/web/. It is active. You can get data using wfs services.
Up to 3 weeks ago I could use QGIS and add wfs layers from my geoserver (http://webservice.dftrans.df.gov.br:8080/geoserver/web/). Now, I can't get layers from it.
I connect to WFS server, I can see the layers, but when I try to load a layer, after timeout, I read this WFS message in my QGIS display:
Network panel:
WARNING    Network request http://webservice.dftrans.df.gov.br:8080/geoserver/wfs?SERVICE=WFS&REQUEST=DescribeFeatureType&VERSION=2.0.0&TYPENAMES=DFTrans:Camada_de_Pontos_de_Rede&TYPENAME=DFTrans:Camada_de_Pontos_de_Rede timed out
WFS panel:
2021-03-03T15:29:53     WARNING    DescribeFeatureType network request failed for url pagingEnabled='true' preferCoordinatesForWfsT11='false' restrictToRequestBBOX='1' srsname='EPSG:31983' typename='DFTrans:Camada_de_Pontos_de_Rede' url='http://webservice.dftrans.df.gov.br:8080/geoserver/wfs' url='http://webservice.dftrans.df.gov.br:8080/geoserver/wfs?request=GetCapabilities' version='auto': Download of feature type failed: Operation canceled
2021-03-03T15:29:53     WARNING    Download of feature type failed: Operation canceled
Messages panel:
2021-03-03T15:29:53     CRITICAL    Layer is not valid : The layer pagingEnabled='true' preferCoordinatesForWfsT11='false' restrictToRequestBBOX='1' srsname='EPSG:31983' typename='DFTrans:Camada_de_Pontos_de_Rede' url='http://webservice.dftrans.df.gov.br:8080/geoserver/wfs' url='http://webservice.dftrans.df.gov.br:8080/geoserver/wfs?request=GetCapabilities' version='auto' is not a valid layer and can not be added to the map. Reason:
At geoserver log file, 'wraper', I can read:
INFO   | jvm 1    | 2021/03/03 15:29:04 | 03 Mar 15:29:04 INFO [geoserver.wfs] -
INFO   | jvm 1    | 2021/03/03 15:29:04 | Request: getServiceInfo
INFO   | jvm 1    | 2021/03/03 15:29:04 | 03 Mar 15:29:04 INFO [geoserver.wfs] -
INFO   | jvm 1    | 2021/03/03 15:29:04 | Request: getCapabilities
INFO   | jvm 1    | 2021/03/03 15:29:04 |     acceptVersions:
INFO   | jvm 1    | 2021/03/03 15:29:04 |         version[0] = 2.0.0version[0] = 1.1.0version[0] = 1.0.0
INFO   | jvm 1    | 2021/03/03 15:29:04 |     baseUrl = http://webservice.dftrans.df.gov.br:8080/geoserver/
INFO   | jvm 1    | 2021/03/03 15:29:04 |     service = WFS
INFO   | jvm 1    | 2021/03/03 15:29:16 | 03 Mar 15:29:15 INFO [geoserver.wfs] -
INFO   | jvm 1    | 2021/03/03 15:29:16 | Request: getServiceInfo
INFO   | jvm 1    | 2021/03/03 15:29:16 | 03 Mar 15:29:15 INFO [geoserver.wfs] -
INFO   | jvm 1    | 2021/03/03 15:29:16 | Request: describeFeatureType
INFO   | jvm 1    | 2021/03/03 15:29:16 |     service = WFS
INFO   | jvm 1    | 2021/03/03 15:29:16 |     version = 2.0.0
INFO   | jvm 1    | 2021/03/03 15:29:16 |     baseUrl = http://webservice.dftrans.df.gov.br:8080/geoserver/
INFO   | jvm 1    | 2021/03/03 15:29:16 |     typeName[0] = {10.233.30.62/Geo_Data_DFTrans}Camada_de_Pontos_de_Rede
INFO   | jvm 1    | 2021/03/03 15:29:16 |     outputFormat = text/xml; subtype=gml/3.2
What I must to do?


texts: Proper name for area defined by coordinate pair and radius

texts:
In defining geography objects within software, I have things such as a City Area, which is an area defined by the bounds of a city, or a County Area, defined by the county boundary.

What is the proper term for an area defined by a center coordinate pair with a radius?

I thought circle might be appropriate, but wasn't sure since that circle is bent along the earth's surface.  Is there a proper name for this?


texts: How does ArcMap decide ArcGIS Online and Basemaps are not reachable?

texts:
I am currently at a remote location and even though I have internet access (i.e. I can browse the web but not very quickly), when I start ArcMap 10.1 the Add Data options for Add Basemap and Add Data from ArcGIS Online are greyed out.

Does anyone know why this might be, and whether there is a way I can force ArcMap to "try harder" to reach ArcGIS Online?


texts:
I have data with columns x,y, address.
And no any information about the projection of x,y points.

How can I convert these points to correct lat/long pairs?

1091969,199515,27 HEDGEWAY CT
1092046,199590,4 SARATOGA CIR
1091986,199663,6 SARATOGA CIR
1091928,199696,8 SARATOGA CIR
1091883,199758,2 LEXINGTON CIR
1091948,199825,4 LEXINGTON CIR
1091917,199913,6 LEXINGTON CIR
1091839,199946,8 LEXINGTON CIR
1092026,198875,20 WARNER AVE
1092013,198934,28 WARNER AVE
1092000,198992,32 WARNER AVE
1091988,199050,WARNER AVE
1091974,199111,48 WARNER AVE
1091953,199210,56 WARNER AVE
1091940,199269,62 WARNER AVE
1087464,233134,7 EDGEWOOD LN

I've tried to use GeoPandas.

Create dataframe as epsg: 2236 crs, because its Florida , Nassau county.
Then convert to epsg:4326 (lat/long CRS)
Then I've checked a few points in Google Maps and tried to adopt just add difference in latitude and longitude as Point(p.x + 6.07406, p.y + 15.83704)

But it seems such a linear approach does not work.

Any advice?

    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.coordinate_x, df.coordinate_y))
    gdf.crs = {'init': 'epsg:2236'}
    gdf['xy_geometry'] = gdf['geometry']

    # convert to latitude/longitude coordinate reference system
    gdf.to_crs({'init': 'epsg:4326'}, inplace=True)
    gdf.rename(columns={'geometry': 'lat_long_geometry'}, inplace=True)

    # correct geometry points with latitude & longitude offsets
    gdf.lat_long_geometry = gdf.lat_long_geometry.apply(lambda p: Point(p.x + 6.07406, p.y + 15.83704))


texts: Displaying dynamic text in Data Driven Pages from joined MS Access database table?

texts:
I'm running ArcMap 10.4 and trying to incorporate text from a joined MS Access 2013 table (through an OLE DB connection) as a dynamic text element in a data driven pages setup. Adding dynamic text from those fields returns "[empty]," even though I can label from those fields and run symbology from those fields. When I export my Access table to Excel and join that instead, dynamic text from those fields works splendidly. I have tried the usual "tablename.fieldname" workarounds and the such to no avail. I've attached a screenshot of the attribute table (with a few fields cropped and/or blurred for confidentiality and size). It's PropertyAddress that I'm attempting to make a dynamic text block for.

The ESRI help specifically says, "this includes fields joined to the index layer," but I'm just not able to make that happen when the fields are joined from Access. Is this a deficiency in ArcMap or have I missed an important memo somewhere?


texts:
Created a local WFS layer with GeoServer.

Loading in ArcGIS Pro seems, X and Y turned and rounded!??

Features are all in Holland with 4,xxxx and 51,xxxxx

So now all the records come in rounded 51,00000 and 4,00000.

QGIS adding same WFS is working well.

Blows my mind.

Someone an idea?


texts: Semi Automated Classification plugin (QGIS) error when trying to create ROI

texts:
When I try to create a new ROI through Semi-Automated Classification plugin (SCP) for QGIS, I keep getting a message saying:

"It appears that SciPy is not correctly installed. Please, see this page for information about SciPy installation  http://semiautomaticclassificationmanual-v3.readthedocs.org/en/latest/installation_ubuntu.html#qgis-download-and-installation"

(the page on this link didn't help me, unfortunately)

Any idea on how to solve this?

Versions:

QGIS 2.8.2

SCP 4.3.3

Mac OS X 10.7.5


texts: How to rearrange FID numeration programmatically?

texts:
I've got a feature class

Now, I delete the selected features programmatically using pFeature.Delete :

Is there any command or function in VBA which rearranges the FID field and makes the numeration to begin from 0 like in the picture below?

Or may be it is possible to delete the features programmatically  without corrupting FID numeration, like the manual deleting do it.

Starting and stopping editing doesn't solve the problem (the code deletes selected features of the first layer):

Dim pMxDoc As IMxDocument
Dim pFeatureLayer As IFeatureLayer
Dim pFeatureSelection As IFeatureSelection
Dim pSelectionSet As ISelectionSet
Dim pFeatureCursor As IFeatureCursor
Dim pFeature As IFeature
Dim pDataSet As IDataset
Dim pFC As IFeatureClass
Dim pUID As New UID
Dim pEditor As IEditor

Set pMxDoc = ThisDocument
Set pFeatureLayer = pMxDoc.FocusMap.Layer(0)
Set pFeatureSelection = pFeatureLayer
Set pFC = pFeatureLayer.FeatureClass
Set pSelectionSet = pFeatureSelection.SelectionSet
Set pDataSet = pFC

pSelectionSet.Search Nothing, False, pFeatureCursor
Set pFeature = pFeatureCursor.NextFeature

pUID = "esriEditor.Editor"
Set pEditor = Application.FindExtensionByCLSID(pUID)
pEditor.StartEditing pDataSet.Workspace
pFeature.Delete
pEditor.StopEditing True

If you execute the Sub below then the FID numbers will be reaaranged! In the end, solution is to start and stop editing in another Sub (not in the same).

Private Sub rearrangeFID()

Dim pMxDoc As IMxDocument
Dim pFeatureLayer As IFeatureLayer
Dim pDataSet As IDataset
Dim pFC As IFeatureClass
Dim pUID As New UID
Dim pEditor As IEditor

Set pMxDoc = ThisDocument
Set pFeatureLayer = pMxDoc.FocusMap.Layer(0)
Set pFC = pFeatureLayer.FeatureClass
Set pDataSet = pFC

pUID = "esriEditor.Editor"
Set pEditor = Application.FindExtensionByCLSID(pUID)
pEditor.StartEditing pDataSet.Workspace
pEditor.StopEditing True

End Sub


texts: Reclass for vector data in ArcGIS Desktop?

texts:
I've got a vector dataset (in fact, OS MasterMap data) and I want to set the value of a certain field based on the value of another field. For example:

Field1        Field2
Building      20
Grass         50
Water         19

In the above, I want to set the value of Field2 based on the value in Field1, that is - all records that have Building in Field1 should have Field2 set to the value of 20.

I know I can do this by doing Select by Attribute and then using the Field Calculator to give the selected records a value, but doing this loads of times is going to get very boring.

Is there a way to give ArcGIS Desktop a sort of table of values like you do in the Raster Reclass tool, and get it to do it automatically?


texts: Translating X, Y coordinates to lat/long

texts:
I have X and Y coordinates from a resized image that has its origins from a Sentinel 2B photo. The original photo is 109.681 by 110.334 km. The re-dimensioned file is 370 by 369. I also know the latitude and longitude coordinates from the upper right and left corner, center and lower right and left corner.
So by multiplying a desired X and Y pixel by 300m(how much 1 pixel equals to) and applying the Euclidean distance formula I can find out how far said pixel is from my reference point(I'm using the upper right corner as reference.)
My first attempt at translating the Cartesian coordinates to lat and long went as following:
def pixel_meter():
    global x_meter
    global y_meter
    x_pixel = int(input('Which X pixel? '))
    y_pixel = int(input('Which Y pixel? '))
    x_meter = (x_pixel * PIXEL_METER) * -1
    y_meter = (y_pixel * PIXEL_METER) * -1

def conversion():
    dLat = y_meter / EARTH_RADIUS

    dLon = x_meter / (EARTH_RADIUS * np.cos(np.pi * LAT_REFERENTIAL / 180))

    lat_final = LAT_REFERENTIAL + dLat * 180 / np.pi
    lon_final = LON_REFERENTIAL + dLon * 180 / np.pi

    print(lat_final, lon_final)

That got me close, but when I checked it on a map online it was still off by some kilometers to what would make sense by judging from the pixels displayed by the resized picture.
Then after searching some more I found about the Vincenty's formula and how more precise it is. With pygeodesy I tried this:
from pygeodesy.ellipsoidalVincenty import LatLon

p = LatLon(-18.088506124573, -44.907758947306) #Reference point
d = p.destination(85854, 224.47) # Distance to my desired X, Y coordinates and bearing from center.
print (d.lon, d.lat)

The distance is calculated using Euclidean's distance.
The result is still off, and I'm certain it is because of the bearing. I don't know, and found no way to know what would be the bearing of my points of interest(more than 4 thousand X,Y points).
Is there any way to calculate the bearing just from the information I have?
That is, latitude and longitude from center, upper corners, lower corners and original covered area. I'm open to any possible methods, really.
Also, if it's not possible to calculate the bearings then is there any other way to make this translation between coordinates more precise?


texts: Creating neat map in basemap using simplified boundaries?

texts:
Is it possible to create a neat map like using basemap (mpl_toolkits.basemap)? I am particularly curious to know if there's a way to "simplify" the boundaries as they are rather edgy in the current version on basemap?
This neat map is borrowed from cartopy http://scitools.org.uk/cartopy/docs/latest/examples/features.html


texts: Plotting points using x,y coordinates from point on map in QGIS?

texts:
In my line of work, fieldworkers measure specific drilling points and give us x,y coordinates from easily identifiable reference points (the corner of a house for example).

Is it possible to plot their reference point on a map (say a satellite image), define an axis and then insert the x,y coordinates they provided (in meters) to find the point they plotted on the map?

I'm using QGIS.


texts:
The field of view analysis in GIS is a polygon file that shows the portions of an area that are visible versus no visible across a terrain starting from a given point. see this Link

How can I use Field of view analysis in ArcGIS? is there any tools to use field of view analysis in ArcGIS or any scripts ?


texts: Which is the simple spatial database for GIS which can is easy to learn?

texts:
I am looking for GIS DB, to work with 800 millions of GIS data.
I want to make my work more simple easy by doing it in spatial DB. I want to work with Geography and not on Geometry.
Please suggest some good and easy to understand the db. I use MAPINFO for mapping.


texts: Does ArcGIS Pro create Blank project always open Project View now?

texts:
Prior to ArcGIS Pro 1.4 my recollection is that when choosing to create a Blank project the Project View did not open - see below for the old behaviour where I have the Contents and Project Panes open but, if closed and the project saved, then they do not open for the next Blank project that is created.  This old behaviour is the same as that seen in the Create a project Help's Tutorial summary video.

At ArcGIS Pro 1.4 creating a new Blank project always has the Project View open.

If this change in behaviour is reproducible between versions, then is it considered to be an enhancement or a minor bug?


texts: Applying filter to features drawn from ArcGIS Server WMS service,?

texts:
I have an ArcGIS 10 Server with various map services running both the REST API and WMS servers. In the REST API I can call the layerDefs argument to filter features used by one of their properties (e.g. area > 50000).

How do I do this with a WMS call?


texts: Reprojecting points in NAD 1983 to LLA (WGS84)?

texts:
I've read documentation about using NAD_1983_To_WGS_1984_1 for the entire North American continent or NAD_1983_To_WGS_1984_5 for the 48 contiguous United States but what I'm confused is how can I use these when my data is in IL plane because 0,0,0 is the corner of IL whereas conventional NAD83/WGS84's 0,0,0 is the equator?

For example, if I have a point 1012111, 1276111 in the IL plane how do I turn it to LLA?

Data source:

PROJCS["NAD_1983_HARN_StatePlane_Illinois_East_FIPS_1201",
    GEOGCS["GCS_North_American_1983_HARN",
        DATUM["NAD83_High_Accuracy_Regional_Network",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",300000],
    PARAMETER["False_Northing",0],
    PARAMETER["Central_Meridian",-88.33333333333333],
    PARAMETER["Scale_Factor",0.999975],
    PARAMETER["Latitude_Of_Origin",36.66666666666666],
    UNIT["Meter",1],
    AUTHORITY["EPSG","102271"]]


texts: How do I prevent Qgis from being detected as "not responding" when running a heavy plugin?

texts:
I use the folowing line to inform the user about the status:

iface.mainWindow().statusBar().showMessage("Status:" + str(i))

The plugin takes about 2 min to run on my data set but windows detects it as "not responding" and stop showing the status updates. For a new user this is not so good since it looks like the program have crashed.

Is there any work around so the user is not left in the dark regarding the status of the plugin?


texts: Using InsertCursor and SearchCursor with Nested Loops

texts:
I have a script that's parsing through text within each record of a table/field, and inserting those values into field/records of another table. Here's an example of the text being parsed (from "Field1"):

Name Multiple Words (0/24.56) Name Multiple Words2* (96/24.56) Name Multiple Words3* (0) Name Multiple Words4* (96/12.58) Name Multiple Words5* (96/78.12) Name Multiple Words6* (0/98.32) Name Multiple Words7* (96/0) Name Multiple Words8* (0) Name Multiple Words9**

Here's the script:

import re, arcpy, sys

# Local variables:
Table2 = "D:\Source_Data_Convert.gdb\Table2"
RAW_Data = "D:\Source_Data_Convert.gdb\RAW_Data_Table2"

#Create Cursors and Insert Rows
insertcursor = arcpy.da.InsertCursor(Table2, ["Other_Stuff", "That_Stuff", "Some_More", "From", "To", "Num_One", "Num_Two", "Nums"])
with arcpy.da.SearchCursor(RAW_Data, ["Field1", "Other_Stuff", "That_Stuff", "Some_More"]) as searchcursor:
########  Start the SearchCursor Loop ###################
        for row in searchcursor:
            try:
                Other_Stuff = row[1]
                That_Stuff = row[2]
                Some_More = row[3]
                listFrom = re.split(r'\*\s*\(.*?\)\s*', row[0])
                print listFrom

                Nums = re.findall(r'\(([^)]+)\)', row[0])
                for match in re.finditer(r'\(([^)]+)\)', row[0]):
                    parts = match.group(1).split('/')
                    print parts
                    First_Num = parts[0]
                try:
                    Second_Num = parts[1]
                except IndexError:
                    #Second_Num = None
                    Second_Num = 0
                    print "First_Num, Second_Num: ", First_Num, Second_Num

                print "Parsing Successful"

############  Start the Insertcursor Loop ############
                for n,Value in enumerate(match): #enumerate is essentially doing a count
                    insertcursor.insertRow((Other_Stuff, That_Stuff, Some_More, listFrom[n], listFrom[n+1], First_Num, Second_Num, Nums[n]))
                    print "Data Inserted"
            except:
                pass
            else:
                break

del insertcursor
del searchcursor
del row

I can get the values to print correctly, but I can't seem to get the (2) values from the inner nested loop ("First_Num" and "Second_Num") to insert correctly with each (1) value ("Seg_Nums", etc.) from the outer loop. When I run the script, the last set of values from the inner loop get populated with each Seg_Num value from the outer loop.  Hopefully my question makes sense. Basically I'm trying to populate (3) fields. One field gets the values from the outer loop, the other two fields get the values (resulting from the split) from the inner loop.


texts: Handling mouse events for tool called by button?

texts:
I have a dockable window with a 2 buttons that call 2 different tools (Draw Marker and Select Features), What i want to happen is when the user clicks on the relevant button on the dockable window and then clicks on the map, a onMouseUp event is called and it runs code i need to run.

Where does the onMouseUp event need to be placed so this can be done?


texts: Creating all possible line segments between all points using QGIS

texts:
I want to connect all points in a point feature with every possible line segment using QGIS.
For example, how a pentagon has five points with ten line segments connecting them all to form a pentagram inside.
Alternatively phrased: How do I make spokes when all the points are hubs?
I don't care about the order, I just want every connection between two points as one line each. Is there a tool or script that can do this in QGIS?
This topic (How to generate line segments between all points) is the closest to what I want but the plugins are no longer supported.


texts: Translating lat/lon coordinate on rotated/scaled maps

texts:
I need to indicate a lat/lon coordinate in an old map like this:
http://www.davidrumsey.com/luna/servlet/detail/RUMSEY~8~1~3017~90020003:Topographical-Map-Of-The-City-and-C

Projected on Google Maps:
http://rumsey.geogarage.com/maps/g2784000.html

I know the bounding box of the map as well as the dimensions of the image in pixels, but I couldn't figure out how to locate a lat/lon within this map.

Upper left is (40.698291,-74.079051), Lower Left (40.659855,-73.979638), Lower Right (40.855232,-73.835582) and Upper Right (40.882919,-73.940263).

Is there a standard formula on how to do locate a lot/lon datum in the historic map?


texts:
I have a PostgreSQL database with "UTF-8" encoding which contains data in Russian language.
Each table of database loaded in QGIS as a PostGIS layer.
I want to show labels of objects on a map.
When I open .qgs file on my machine in QGIS 3.2.2 it works fine.

But when I open it via QGIS Server 3.2.2 on dedicated machine I have empty squares instead of text.

Both machines work under Windows 10 (English)
What the problem could be?
Update1:
I've also installed QGIS 3.2.3 on my dedicated machine. Russian text also doesn't work there. May this problem be in Apache configuration?
Update2:
I've created another win-1251 database. It also doesn't work.


texts: Can any suggest an open source workflow to convert an Esri LiDAR terrain dataset to a LAS file

texts:
I’ve got a very large file geodatabase (1TB) that contains LiDAR data stored as a multiple terrain datasets.

I’d like to spatially filter the data and export the result as a LAS file without using Esri tools (because I don’t have  access to them).

I’ve previously dumped sections as xyz files, the. Used LAStools to covert to LAS - but it’s extremely time consuming, which doesn’t really justify the lastools licence (to be clear - the LAStools components are extremely fast, it’s the filegdb to xyz that’s slow).

Has anyone got a filegdb to PostGIS/PDAL/las workflow that works?


texts:
I am working on ArcGIS Online.I want to add New Labels to seperate my fields.
The first image represents the original Fields.

The second image represents how I want it to be(Seperated with labels).

How can I add Labels to ArcGIS Online to seperate my fields?


texts: Type of event dispatched when the user double-clicks the map

texts:
I am currently working with a drawing tool for the ESRI Flex API, and every time I double-click the mouse a map service will perform a measurement and display the length of the line that I am drawing.

I want to mimic this double-click manually by dispatching a MouseEvent.DOUBLE_CLICK, but the map service is not listening to this default Flex event. I suspect that the API has created a custom MapMouseEvent or something like it that is being dispatched when a user double-clicks the mouse.

Does anyone know which event is being dispatched when I double-click the mouse on the map?


texts: Georeferencing HiRISE images in QGIS 3.6 on Windows

texts:
I have a CTX mosaic for a base layer in my project. However, I am running into trouble overlaying HiRISE images over it. HiRISE final products should be georeferenced, but they just do not fall into place. I tried the manual georeferencing tool in QGIS to fix it but it also doesn’t work.
An example of a HiRISE image is https://hirise-pds.lpl.arizona.edu/PDS/RDR/ESP/ORB_036900_036999/ESP_036925_1985/


texts: Why does TableToTable give ERROR 000732 when published to ArcGIS Server?

texts:
I have a python code that uses ArcPy command "arcpy.TableToTable_conversion" to convert a database table to a dbf file.

arcpy.TableToTable_conversion

The input table is from a SQL Server database through an odc Database Connection in ArcCatalog.

"Database Connections\\myDb_DEV.odc\\dbo.V_Locations"

When run straight from ArcCatalog, the code worked perfectly. However, after it was published as a Geoprocessing service on ArcGIS Server 10, it failed to run.

Error
Error Executing Task. Error executing tool.: : Failed to execute. Parameters are not valid. ERROR 000732: Input Rows: Dataset "Database Connections\myDb_DEV.odc\dbo.V_Locations" does not exist or is not supported Failed to execute (TableToTable). Failed to execute (ExportToBOMAWS).

Code: 500

Could anybody help with any idea why and how to make it work as geoprocessing service?


texts:
I have a webservice where I have one solution to obfuscate location called Points of interest which are predefined points scattered around an area.

When ever a user send their coordinates from their cellphone I calculate the Euclidean distance from their current location to the closest point of interest resulting in their location data is associated with the predefined point of interest.

This is fine, but I need one more way to obfuscate location.

What other options are there to obfuscate location?

Can I do something so simple as trimming the LatLng to two decimals?

How about adding noise to the coordinate data?

What type of noise can I use, will gaussian perturbation work?

I read about perlin noise https://en.wikipedia.org/wiki/Perlin_noise


texts: Gdal: How to get a pixel's value from raster image?

texts:
Given a topographic raster image.

Given a pixel's coordinate in pixel/lines or geocoordinate.

How to get it's pixel value ?