texts:
This has been bugging me for a while. I can easily call a form from a button using ShowDialog() if I use standard vb.net windows form. However when I use the ESRI addin dockable window which uses hooks in the form I can't seem to open it.

It requires a 'hook as object' when I try to call it from a button. ESRI forms do not use the ShowDialog() method and require you to use Show(). I have searched and look everywhere but can't see how to call the form that uses hooks. Do I need to create a new hook object in the button class?

When I click the button nothing happens and it turns grey. The form is not hidden behind ArcMap or anything.

Here is what I use to call the form:

Protected Overrides Sub OnClick()
    Dim oform As New window(Me.Hook)
    oform.Show()

End Sub

Here is the ESRI form

Public Class window

Public Sub New(ByVal hook As Object)

     InitializeComponent()

    Me.Hook = hook
End Sub

Private m_hook As Object

Public Property Hook() As Object
 Get
  Return m_hook
End Get
Set(ByVal value As Object)
  m_hook = value
End Set
End Property

Public Class AddinImpl
Inherits ESRI.ArcGIS.Desktop.AddIns.DockableWindow

Private m_windowUI As window

Protected Overrides Function OnCreateChild() As System.IntPtr
  m_windowUI = New window(Me.Hook)
  Return m_windowUI.Handle
End Function

Protected Overrides Sub Dispose(ByVal Param As Boolean)
  If m_windowUI IsNot Nothing Then
    m_windowUI.Dispose(Param)
  End If

  MyBase.Dispose(Param)
End Sub

  End Class

End Class

UPDATE:

Now I tried

Dim dockWindow As ESRI.ArcGIS.Framework.IDockableWindow

    Dim dockWinID As UID = New UIDClass()
    dockWinID.Value = My.ThisAddIn.IDs.window
    dockWindow = My.ArcMap.DockableWindowManager.GetDockableWindow(dockWinID)

    dockWindow.Show(True)

And still does not work.


texts: Adding band to GDAL VRT dataset and converting to TIFF using GDAL?

texts:
I'm trying to add a band to a TIFF dataset by converting it to a VRT format and using the AddBand() function, then convert it back to a TIFF to access the band and write the data. Here is an example of my code:
const char *formatVRT = "VRT";
const char *formatTIFF = "GTiff";
char *filenameTIFF;
char *filenameVRT;
GDALDriver pDriverTIFF, pDriverVRT;
GDALDataset pTiffDS, pVrtDS;

GDALALLRegister();

pDriverTIFF = GetGDALDriverManager()->GetDriverByName(formatTIFF);
pDriverVRT = GetGDALDriverManager()->GetDriverByName(formatVRT);

pTiffDS = pDriverTIFF->Create(filenameTIFF, cols, rows, bands, GDT_Float32, NULL);
pVrtDS = pDriverVRt->CreateCopy(filenameVRT, pTiffDS, FALSE, NULL, NULL, NULL);

pVrtDS->AddBand(GDT_Float32, NULL);
pTiffDS = pDriverTIFF->CreateCopy(filenameTIFF, pVrtDS, FALSE, NULL, NULL, NULL);

I get the following errors:

ERROR 1: Deleting (filenameTIFF) failed:
Permission denied
ERROR 1: TIFFFillStrip:Read error at scanline 25; got 0 bytes, expected 24840


texts: Render order of vector features in OpenLayers 3

texts:
In a web map application I have a vector layer from a GeoJSON source of country borders. The layer loads perfectly and interactions such as ol.interaction.select works fine as well.

I want to style individual features from the layer when the map loads and this is where I have some problems. The features in the GeoJSON file is ordered according to the id property, in ascending order. Since the id of the features to be styled are retrieved from a database, similar order (or id) of the features in the db and in OpenLayers is crucial.

In OpenLayers I now try to retrieve features from the GeoJSON source assuming that the order of the features is the same as in the file, but the ordering of the features seems completely random. Since the features are not assigned id in OpenLayers I loop through the array of features from the source and assign them ids just to be able to compare the ordering between file and OpenLayers.

 map.beforeRender(function () {
     var features = countrySource.getFeatures();

     //Assigning id to features
     for (i = 0; features.length > i; i++) {
          features[i].setId(i);
        }
    });

map.on('singleclick',function(evt){

    //alerts the id of the clicked country, id from file
    var fileId = select.getFeatures().item(0).get('id')
    alert(fileId);

    //alerts the assigned id
    var openLayerId= countrySource.getFeatures()[fileId].getId();
    alert(openLayerId);

});

This gives me completely different id from the two alerts. It seems to me that the order of the features in the resulting array from source.getFeatures() method is random compared to the file.

Here is a line from the GeoJSON file:

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

"features": [
{ "type": "Feature", "properties": { "X": 66.0867, "Y": 33.8564, "id": 0, "country":
"Afghanistan" }, "geometry": { "type": "Polygon", "coordinates": [ [
[ 61.210817091725744, 35.650072333309225 ], [ 62.230651483005886 ......

According to the API there is a 'renderOrder' option when creating vector layers, http://openlayers.org/en/v3.1.1/apidoc/ol.layer.Vector.html?unstable=true

Is this something I can use to control the ordering of the features when loading them into my vector layer? Or is there another way around my problem of styling the right features based on id?


texts: choosing processing extent from a polygon in a feature

texts:
I am currently trying to convert a list of polygons to individual rasters. Using Arcgis Desktop's 'convert polygon to raster', this is rather simple. But since the polygons are rather spread out, each time i create a raster (for 1 polygon), it uses the entire extent of the feature class, which means the resulting raster is unnecessarily big (with many nodata cells).

I was therefore wondering, if its possible to choose the specific polygon's extent, instead of choosing all polygons' extents? I have tried extracting each separate polygon beforehand, but it seems a bit slow.

I will try to implement the entire analysis (which is a bit longer) in arcpy afterwards, so if the solution is easier to find in there, that might help aswell.

edit: I think i might i have been vague in my description. I mean i have a feature class containing multiple polygons. is there a way to select the extent of one of these polygons? right now i can only select layers (i.e. the feature class and all of its polygons), not separate polygons in the feature class.


texts: Rasterise a shapefile with gdal into a new raster file with multiple bands

texts:
I want to convert a shapefile (polygons) which contains also 3 attributes (CO,NOX,SOx) with float values.

I used the following command

gdal_rasterize -a CO -a NOx -a SOx -ot Float32 -ts 1000 1000 -l emissions_shp /home/emissions_shp.shp /home/emissions.tif

and the result is a raster file with one band only. It seems that it gets only the last attribute (Sox).


texts: Is there a way to allow user to change just the beginning of a path name in a toolbox for several parameters at once?

texts:
I am altering someone else's Python code to create a GIS tool that everyone can use.  Currently, they have about 20-30 feature classes with paths to them in their code. The problem is that the beginning of those paths will be different depending on who is using them and where they paste the main folder with the features in it.
For example:
They had
arcpy.MakeFeatureLayer_management(r"V:\gdrs\data\fens\fgdb\biota.gdb\fens", "nearFeature")

I pulled out the path and made the path a parameter
Fens = arcpy.GetParameterAsText(5)

and switched their code to
arcpy.MakeFeatureLayer_management(Fens, "nearFeature")

Then I make it a feature class parameter in the toolbox.
However, as I'm changing this I realized there are about 20-30  features so this isn't practical to have the user browse to the file location of 30 features each time they use the tool.
It would be nice if they could just change the beginning of the path as a parameter. It would also be nice if they didn't have to change the path again after they've used to tool once.
Note: gdrs is the main folder and everything should be the same after that.
I'm new to Python, but I am assuming there is an easy way to get around this issue.


texts: How do I update a polygon with values from a line in GRASS

texts:
I have a line map generated from r.stream.extract.  I also have a polygon map generated from r.stream.basins, using the raster stream layer generated from r.stream.extract.
I want to update the polygon layer with the category and one of the attributes of each line.
Each polygon contains one and only one line, which touches the boundary at either one or two points.
In PostGIS, I'd do something like the following:
UPDATE poly_table a
SET line_cat = (
    SELECT line_cat FROM line_table b
    WHERE st_contains(a.geom, b.geom)
);


texts: How to See The Attribute Result Of Geometric Network Trace

texts: Running ArcMap 10.4.1 I can trace a geometric Network and see the Graphical returns of Junctions and Edges inside the network but can you please let me know How/ Where I can see the Attributes return from the Trace?


texts: What is difference between ArcGIS feature layer and layer?

texts:
In the ArcGIS help documentation, they tend to use the terms "layer" and "feature layer" interchangeably (see for example "What is a layer?" or Make Feature Layer documentation).

Is "layer" simply an abbreviation for "feature layer"? Is there any instance where a layer would not be a feature layer, or vice versa?


texts: Is there an advanced way to label in ArcMap so that the label will not be dropped?


texts: Changing URI source of WMS Raster Layer in place PyQGIS

texts:
In my Qgis Plugin, I manage a set of WMS layers that are secured by Api key. As I develop this plugin for differents clients who have different Api Key, I need to set the right Key when a specific client connect to the plugin.
As I can change the DataSource of a vector Layer without reloading it and keeping it in place in the project with the setDataSource() method of QgsVectorLayer, I wonder if I can do the same with a QgsRasterLayer.
I can't find a method in Api docs.

For now my workaround is to remove all my WMS Layer and reload it from scratch with the right Api Key. For me it's not optimal because it doesn't keep the layer tree clean... (layer order, groups,...) and it takes time...

Here is the code that I actually use:

layers = [layer for layer in QgsMapLayerRegistry.instance().mapLayers().values() if layer.name() in all_flux_provider]
# all_flux_provider is the list of WMS layer that I manage
        for layer in layers:
            # remove all managed flux
            QgsMapLayerRegistry.instance().removeMapLayer(layer.id())
        #Add client layers from scratch
        for flux in fluxs: # fluxs is a list of dict that represent the WMS layer of the client (it comes from a database). structured like that:
                           #{key:apikey of the client, provider:name of the flux provider, url:uri to load the flux}
                layer = None
                if '{insee}' in flux['url']: # {insee} is a placeholder in the url of the flux
                    for insee in insees:
                        uri = flux['url'].strip().format(insee=insee)
                        layer = QgsRasterLayer(uri, flux['provider']+'_'+str(insee), 'wms')
                        QgsMapLayerRegistry.instance().addMapLayer(layer)
                elif '{client_key}' in flux['url']: #{client_key} is a placeholder in the url of the flux
                    uri = flux['url'].strip().format(client_key=flux['key'])
                    layer = QgsRasterLayer(uri, flux['provider'], 'wms')
                    QgsMapLayerRegistry.instance().addMapLayer(layer)

Does anyone know a better way or a good trick to achieve this?


texts: Raster to vector conversion using GDAL Python?

texts:
I have used the following code and it gives me the output shapefile which is of like 4.36 MB but when I open it in QGIS to view the shapefile it shows nothing:

from osgeo import gdal, ogr
import sys
import os

gdal.UseExceptions()
os.chdir("D:/Nikka")
fileName = "D:/subset_raster.tif"
src_ds = gdal.Open(fileName)
if src_ds is None:
    print('Unable to open %s' % src_fileName)
    sys.exit(1)
srcband = src_ds.GetRasterBand(1)
dst_layername = "PolyFtr"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None)
newField = ogr.FieldDefn('Area', ogr.OFTInteger)
dst_layer.CreateField(newField)
gdal.Polygonize(srcband, None, dst_layer, 0, [],
callback=None )

for reference GDAL polygonize in python creating blank polygon?


texts: How to fill the area bounded by lines in QGIS?

texts:
There are two types of area I want to fill color:

The areas bounded by a line, such as a building
The area bounded by different lines, such as the river

I don't know how with QGIS. I only know how to fill colors in polygons but bot lines .May you provide some suggestions?

Thanks!

The image here shows that the building and rivers are composed of lines:


texts: Writing Python script for Altering Multiple Fields

texts:
I'm a novice programmer and I would like to ask for something which will facilitate my work. I need to rename all the field names in some shapefiles I am working on and since renaming is an arduous task, I need to automate using some scripting. Renaming a shapefile field name is not possible so I used geodatabase. Now, the 'Alter Field' tool is limited to only one field name at a time.

Can somebody help me with the script I need?

Below is the example.

Example:

arcpy.AlterField_management(r'D:\Workspace\11-30-2014\Approved GDB\Nisma_GIS DATA\2NISMA_WETCO.gdb\NISMA_WETCO\CHAMBER2', 'LOCATION_I', 'LOCATION_ID', 'LOCATION_ID')

How can I include these: ('X_COORD', 'X_COORDINATE', 'X_COORDINATE', 'Y_COORD' , 'Y_COORDINATE', 'Y_COORDINATE'...etc)

I tried just putting all in the same code but I got this error:

  TypeError: AlterField() takes at most 4 arguments

I'm using ArcGIS 10.2.2, Standard (ArcEditor) License.


texts: (Resolved) How to export annotations in QGIS?

texts:
Is there a way to access annotations, save or export them as csv or xls?

Where they are stored? I can´t find that data anywhere.


texts: ArcGIS Desktop in VMWare not connecting to USB?

texts:
I'm using ArcGIS Desktop 10.5 in a VMWare environment (university provided) with VMWare set to open USB when inserted.  I have inserted my thumb drive and have successfully worked within the folders, but upon trying to begin my mapping, data will not add.  I tried to set the path to the folder in ArcCatalog, but I receive a message saying the path is an "invalid or non-existent directory" despite being able to see these pathways in the connecting window.

I've connected and opened the files in the VMWare, it seems that this is an ArcGIS problem, not something about my VMWare settings nor the host computer settings, which, by the way, is a mac.


texts: Clipping Vector File based on certain Raster Values using ArcGIS Desktop?

texts:
I need to clip a vector file based on certain values of a raster (Slope > 2%). I can´t find a direct way to do it.

Trying to vectorize the raster first didn´t work.

Is there an ArcGIS tool I could use that I haven´t found yet?


texts: Modify MapCruncher sample to make it generate tiles with GCS 1984 spatial reference

texts:
I want the MapCooker sample to generate tiles with a projected spatial reference of GCS WGS 1984.  I currently initialize it with a path to a layer with the NAD_1983_HARN_Adj_MN_Isanti_Feet coordinate system.   The sample is here:

http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#/Sample_Multithreaded_MapCruncher/000100000mmm000000/

I already modified the sample to produce .pngs, by changing the StorageFormat to

cacheStorageInfo.StorageFormat = esriMapCacheStorageFormat.esriMapCacheStorageModeExploded;

Now I want to control the projected .pngs spatial reference.  Any idea how to do this?  I've been messing around with it for 2 hours with no luck.

I'm using ArcEngine 10 with C# and VS2010


texts: Update to Bonn 3.2 QGIS seems to have had a effect on some of my code

texts:
I asked a question here:Adding to attributes based on shared line in QGIS? but it seems that in Bonn 3.2 QGIS this does not work.  Any ideas what has changed for it not to work, other snippets of code seems not to work as well since the update?  For example, this is the code from the post above:

SELECT t0.id1 AS id, Group_Concat(t0.id2, ',') AS Adjoin_Po
FROM(
  SELECT t1.fid AS id1, t2.fid AS id2
  FROM Grid AS t1, Grid AS t2
  WHERE t1.fid != t2.fid
  AND st_length(st_intersection(t1.geometry, t2.geometry)) > 0
) AS t0
GROUP BY  t0.id1


texts: Create new geodatabase in modelbuilder and place queried data in new database

texts:
I have built a model that will take an existing database and query all feature classes in that dataset using a SQL Expression and I then want all the queried data from each feature class to be exported by feature class to a newly created database.  I know I can manually create a new geodatabase and then set the output to that the geodatabase in the "Select" tool but I want the model to create the database automatically.  The model here shows how I have tried to do it but by doing it the way illustrated here it creates a new GDB for each iteration, exports the feature to the new GDB and then moves on the next iteration, overwriting the previous feature class.

How can I create the geodatabase first (once and only once), then run the iteration and place each each iterated feature class into the newly created database?

Any advice?


texts: Runtime split our the line based on moving nearest vertex in pgrouting?

texts:
I attached my screenshot using pgrouting.
Is it possible to runtime split our the line based on moving nearest vertex in pgrouting?
Is there any easiest function or solution over there?


texts: Lines sometimes disappear when editing Postgis layer in QGIS 2.16

texts:
I have a GIS tech editing a watercourse layer that is located in a Postgres/PostGIS table. She is using QGIS 2.16 on Windows 7 and I am running PG 9.5. All editing, which includes modifying existing watercourses and adding new one, saves fine in the DB. There are no permission issues, no errors, no problems whatsoever.

However, every now and then she adds a new line, presses save, turns off the editing session, and the line disappears. I checked the PG logs and the QGIS error console and no errors (relating to this issue). The only repeated PG log error relates to permission issues relating to the topology schema – which is probably another issue altogether. Bottom line: at the time of these events, there are no logged errors or QGIS errors.

Does anyone know how to prevent such data loss when creating new rows?


texts: Add points on zoom and remove points when zooming out

texts:
Problem

I would like to add points/shelterIcons when zooming to level 12 or more and remove them again when zooming out. If all the points were defined as a layer I could add/remove them with the function down below. But the points are in GeoJSON files. I don't know how to "handpick" them and put them in a layer.

Code

//=======================
//Define map and map size
//=======================

var map = L.map('map', {
    center: [55.676098, 12.568337],
    zoom: 7,
    minZoom: 2,
    maxZoom: 18,
    zoomControl: false
});

L.control.zoom({
     position:'bottomright'
}).addTo(map);

L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '&copy; OpenStreetMap contributors'
}).addTo(map);

//============
//Shelter icon
//============

var shelterIcon = L.icon({
    iconUrl: 'http://www.clker.com/cliparts/f/e/6/c/120743200840710503sleeping shelter white.svg.thumb.png',
    iconSize:     [20, 20], // size of the icon
    shadowSize:   [50, 64], // size of the shadow
    iconAnchor:   [0, 0], // point of the icon which will correspond to marker's location
    shadowAnchor: [4, 62],  // the same for the shadow
    popupAnchor:  [0, 0] // point from which the popup should open relative to the iconAnchor
});

    //CENTERS VIEW WHEN USER CLICKS ON MARKER
function clickZoom(e) {
    map.setView(e.target.getLatLng(),map.getZoom());
};

//=============================
//Add GEOjson with feature info
//=============================

var RaslangenHalenRundt = L.geoJSON(RaslangenHalenRundt, {
      pointToLayer: function (feature, latlng) {
    return L.marker(latlng, {icon: shelterIcon}).on('click', clickZoom);
  }, onEachFeature: onEachFeature
}).addTo(map);

var Kalmarsundsleden = L.geoJSON(Kalmarsundsleden, {
      pointToLayer: function (feature, latlng) {
    return L.marker(latlng, {icon: shelterIcon}).on('click', clickZoom);
  }, onEachFeature: onEachFeature
}).addTo(map);

function onEachFeature(feature, layer) {
  var layerType = layer.feature.geometry.type;
  if (layerType == 'Point'){
    var PointPopupContent = feature.properties.name;
    layer.bindPopup(PointPopupContent);
  } else if (layerType == 'LineString') {
    var LinestringPopupContent = feature.properties.lengthdistance;
    layer.bindPopup("Route length: " + LinestringPopupContent + " km");
  } else {
    console.log('lengthdistance not present');
  }
}

    //WHEN USERS ZOOMS TO MORE THAN 12 THE SHELTERS APPEAR
map.on('zoomend', function() {
    if (map.getZoom() <12){
            map.removeLayer(shelterIcon);
    }
    else {
        //shelterMarkers.addTo(map);
        map.addLayer(shelterIcon);
        }
});


texts: How do I make a heat or Choropleth map that uses polygons in QGis?

texts:
I've been playing around with the heatmap plugin for QGIS and I am wondering if it is possible to make heatmap with using a similar plugin/tool that uses polygons as well as points?

I would like to get a representation of the density of the data in much the same way that the heatmap plugin does with points.

The following image is a rough sketch of my desired result:


texts: Plot points by longitude and latitude on a map created by centroids

texts:
I'm searching for a way to plot some points for which I have the longitude and latitude (using Google Maps API) on a map created with centroids. I succeeded to do it by using functions such as map_data() and functions from the package ggplot2.

Results (just to show it really works) :

However I don't like the french departments datas given by map_data(), it's quite ambiguous but the nice thing is that their coordinates are given by longitude and latitude. (french departments are similar to counties)

What I would like to do is to use a shapefile from a site such as GEOFLA (here) that provides me datas in projection LAMBERT-93 to plot my map. However the shapefile only contains the coordinates of deparments centroids so, since I have, in one hand, coordinates given by longitude and latitude and, in the other hand, coordinates given by centroids, so I can't plot the points on my map since the coordinates are not in the same "unit".

That's what I get with summary() on the shapefile :

Object of class SpatialPolygonsDataFrame
Coordinates:
        min     max
x   99217.1 1242417
y 6049646.3 7110480
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m
+no_defs]

Summary of map_data() (I show only the coordinates) :

      long               lat
 Min.   :-5.14209   Min.   :41.37
 1st Qu.: 0.05213   1st Qu.:44.41
 Median : 2.45500   Median :46.79
 Mean   : 2.45851   Mean   :46.47
 3rd Qu.: 4.97283   3rd Qu.:48.44
 Max.   : 9.56267   Max.   :51.10

This shows that the coordinates used are different.

I tried many things such as using proj4string(), sptransform() for example to get the same units but I didn't get what I wanted.

So how should I do to get the same units for my two differents types of coordinates ? Is there a way to get the longitude and latitude from the centroids of the shapefile or should I use a function (if it existed) that can transform my points with longitude and latitude into centroids ?

This is my first time using R for maps and polygons, projections etc. also so don't be surprised if I said some real stupid things.

Exemple of points I want to plot :

  centerid latitude longitude                                       adresse_complete
1   0121H1 46.22368  5.210819               900 Route de Paris, 01440 Viriat, France
2   0121H2 46.22368  5.210819               900 Route de Paris, 01440 Viriat, France
3   0218H3 49.37008  3.335127 46 Avenue du Général de Gaulle, 02200 Soissons, France
4   0303H1 46.34292  2.608996       18 Avenue du 8 Mai 1945, 03100 Montluçon, France
5   0303H2 46.34292  2.608996       18 Avenue du 8 Mai 1945, 03100 Montluçon, France

dput() from above :

structure(list(centerid = c("0121H1", "0121H2", "0218H3", "0303H1",
"0303H2"), latitude = c(46.2236804, 46.2236804, 49.3700842, 46.3429172,
46.3429172), longitude = c(5.2108193, 5.2108193, 3.335127, 2.6089958,
2.6089958), adresse_complete = c("900 Route de Paris, 01440 Viriat, France",
"900 Route de Paris, 01440 Viriat, France", "46 Avenue du Général de Gaulle, 02200 Soissons, France",
"18 Avenue du 8 Mai 1945, 03100 Montluçon, France", "18 Avenue du 8 Mai 1945, 03100 Montluçon, France"
)), .Names = c("centerid", "latitude", "longitude", "adresse_complete"
), row.names = c(NA, 5L), class = "data.frame")


texts: How to join point attribute to polyline layer?

texts:
I have a polyline layer in rectangular shapes representing boundaries. These polylines are adjacent to one another. I also have a point layer with the name of the boundary centered within the polyline. I want to join the 'boundary name' attribute field in the point layer to the boundary polyline the point is within.
My method to accomplish this was to polygonize the polyline layer and join the point attributes using the 'within' criteria. However, when attempting to polygonise the polyline layer, the FID of the polyline layer is not retained and the new polygons are assigned new FID's which I was hoping to retain so I could join the data to the original polylines. How would I accomplish joining a text field in a point layer to a polyline layer?
qgis version 3.10.14


texts: With CTM how can I get the latitude and longitude bounding box of a geoPDF?

texts:
Embedded in a geoPDF is info on the CTM - Coordinate Transformation Matrix.

I am unfamiliar with this CTM.

I know how to use python-numpy.

How (or is it even possible) for me to use the CTM information to get the latitude and longitude bounding box of the image?

Is more information need such as the pixel bounding box of the image?

Is there a simple matrix calculation some one can explain to me that I can do in numpy?


texts:
Given an EPSG projection (say, this Alabama one: http://spatialreference.org/ref/epsg/26729/)

How can you take the given WGS84 projection bounds in such a way that you can use them in a D3.js projection.

For example, how would you know what projection, degree of rotation or bounding box to use to show the map?


texts: Seeking geology fill pattern libraries for QGIS?

texts:
I'm making a geological map using QGIS but I don't know how to draw a limestone lithology in QGIS.

To draw other lithologies I combine lines and dots in the best possible way but that is very limited.

Is there a fill pattern library that I can use to do this?


texts: Display multiple attributes for one point

texts:
I am using ArcGIS 10.  How can I display multiple attributes for a single point.  For example, I would like to display (graphically if possible) the amenities provided at various recreational areas... such as my point location "Graham's Landing" and have a symbol for camping, fishing, boat ramp, bathrooms and etc that reference the data in the attribute table.  I would like to do this for each point (rec facility) that I have.

Thanks and I am open to other ideas as well!


texts: How to fix label contour lines layer displays on top of all layers in QGIS

texts:
How can I fix the situation, when the label contour lines  displays (on top of all layers) on top of this layer?

QGIS 2.15.0-Master


texts:
I am trying to overlay the layers that are rendered using MapServer over Google maps.  I am pretty much new to this whole thing. So I tried to go through this example here http://mapserver.org/output/tile_mode.html.

Below is the code which I am using. Could someone let me know where I could be wrong or if there is any other way to configure this.

 <!DOCTYPE html
   PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

 <meta http-equiv="content-type" content="text/html; charset=utf-8"/>
 <title>Google/MapServer Tile Example</title>


texts: QGIS: prepare data for timeline animation

texts:
I have a data-spreadsheet containing a value-attribute that changes over time for a number of administrative units. I have a polygon layer as well with features for each of the administrative units. The spreadsheet containing the date looks like this:

I now want to use the built-in temporal control in QGIS 3.16 to visualise the change of the attribute-value for each administrative unit over time. I imported the spreadsheet as attribute data. Now the problem starts: how to join the attributes to the geometry? Each geometry (administrative unit) should contain for every time-stamp a separate value.
Having the Admin. units in the spreadsheet in the top row, I can't make a table-join with the polygon-layer. So I tried to first convert the spreadsheet in this way:

No I can joint these attributes using the Ad.Unit-field to the geometries of the administrative units. That works, however, I only get the first data-set per ad.unit (day 1). All the others are left out. But to be able to animate, I should have a value for each day (timestamp-field).
What I would like to achieve is a time-animated circle-marker whose size depend on the value-field for each day. How to achieve this?


texts: How to extract elevation values at line start and end points from a DEM?

texts:
I am new to the whole gis-experience and need your help working with QGIS.

I have:

DEM raster with heights (3d)
line-shapefile which contains the rivers (2d)

I need:
heights for "start" and "end" in lines (to get them 3d) then compare that heights with "flow direction" of the lines

qgis-help doesn't really help me out. Tried searching google, nothing useful found. I am pretty sure that it may be a really simple "standard-tool" that can do what I want. However I haven't found it yet.

Edit:

well, i am looking for a solution, not important wich way. dont have any experience with python.

i was working on that problem and think i am almost ready, but there seems to be a problem:

i extracted the nodes of the lines, they got the heights from the dem-grid. i even managed to add the values of the points back to the lines, but when i check the values, the "minimal height" (or maximal) dont always have the suitable x- and y-coordinates of the point.

example:
there is a line (columns: name, id, start_x, start_y, end_x, end_y). i imported the point-values (x, y, z coordinates) into the line. got 6 new columns in the table: min_x, min_y, min_z, max_x, max_y, max_z

now the min and max coordinates are mixed up and dont compare to those from the points.

my thoughts now are:

i delete the min_x, min_y, max_x, max_y
rename "max_z" and "min_z" to "up" and "down"
somehow i need to get the up_x, up_y and down_x, down_y from the points. maybe there is a way to do that in the field calculator like: (for up_x in table "lines") value = "x" FROM "points" WHERE "up.lines" = "z.points"

dont know if that can work in qgis or how the code should look like. maybe there is a way to manage this in python, but as i said, i dont have a clue about python.


texts: I have a shp file with some rectangles and I want to extract statistics within these areas from a raster image. I prefer ArcMap 9.3.


texts:
I am trying to calculate field with the formula. I tried arcpy.calculate field management but still having difficulty.

(10/10.4*(100*(math.exp(0.000871* !A!)-1+math.exp(0.000537*
!B! )-1+math.exp(0.000487* !C! )-1)))

I used the below code.

with arcpy.da.UpdateCursor(fc, (VALUE, A, B, C)) as cursor:
   for row in cursor:
     row[0] = (10/10.4*(100*(math.exp(0.000871* !A!)-1+math.exp(0.000537*
!B! )-1+math.exp(0.000487* !C! )-1)))
     cursor.updateRow()


texts:
I've created a role in GeoServer called WFS_ACCESS. I've assigned it to a user called wfs_viewser and have the following rules setup:

Though when i hit my wfs service and login with wfs_viewer credentials i get an access denied error. As soon as i grant wfs_viewer "ROLE_ADMINISTRATOR" it works fine. What is wrong with my rules?

  Build Information Version
  2.1.0 Subversion Revision 15809 Build Date 11-May-2011 20:09 GeoTools
  Version
  2.7.1 (rev -1)


texts: Issues with osgeo - qgis.core and arcpy modules

texts:
Is it possible to have the python shell from either qgis or argis to load (import) both qgis.core and arcpy?

I am developing a metadata creation tool and have completed the core functionality using arcpy for formats that arcgis reads but want it to handle .tab files as well without having to convert the files to shp and then analyse them as in this post How to read MapInfo .tab files in ArcPy?

===
Thanks to RK I got some of this working (see code below)

BUT I can't figure out how to get the projection and other important info...http://www.gis.usu.edu/~chrisg/python/2009/docs/ogr.pdf includes a GetSpatialReference, GetGeometryType( self ) and GetGeometryName( self ):but how do you use them?

import osgeo, os
from osgeo import ogr
os.chdir(r'L:\\Data_Admin\\QA\\Metadata_python_toolset\\test\\Local_Govt\\TRC\\TAB')
os.getcwd()

driver = ogr.GetDriverByName('MapInfo File')
fn='trc_boundary.tab'
dataSource=driver.Open(fn,0)
layer = dataSource.GetLayer(0)
feature = layer.GetFeature(0)
extent=layer.GetExtent()
numFeatures=layer.GetFeatureCount()
#id = feature.GetFieldAsString('id') # this didn't work
geometry = feature.GetGeometryRef()
x=geometry.GetX()

#delete datasets from memory
#feature.Destroy()
#dataSource.Destroy()

And is there a tool similar to arcpy.ListFeatureClasses() if not I guess we can use glob as anyway for ogr the driver has to be specified.


texts: Perform Calculations in Cartesian Frame (NED/ECEF) Vs. Geodetic Frame

texts:
I'm recording robot's path (using INS) in Cartesian frame coordinates (x,y,z) relative to origin point. So essentially I've created NED coordinate system.

Origin point has exact (in advance measured) geodetic location (LLA: Latitude, Longitude, Altitude).

The calculation I apply for each NED (x,y,z) robot's location to get geodetic position:

Represent the origin point in ECEF.
Because relationship between the origin point to robot's location (x,y,z) is in Cartesian coordinates (NED), and from previous step I have origin point in ECEF, which is also Cartesian. Using well known geometry I can represent robot's location - NED (x,y,z) in ECEF.
Using well known equations I can convert that ECEF location to LLA (geodetic coordinates), for example in WGS84 standard.

Moreover - I can perform spatial calculations in my Cartesian frame coordinates using simple geometry, instead of getting used to complex geodetic calculations.

Question- It sounds too good to be true that one can just using simple geometry and having closed form ECEF to LLA equations to transform into geodetic locations. Assuming perfect transformations (with only numeric errors) and no "Great Circle" issues, where this mode of operation will not be good? For example, if in Cartesian frame the robot moved 100m in X direction, If I transform start & end locations from Cartesian to geodetic, it will still keep 100m distance in WGS84 system?


texts: Solar radiation modelling on vertical surfaces

texts: Are there any methods to measure solar radiation on vertical building facades using Esri's Solar Analyst?


texts:
I compiled the last ledaps from source, but when I try to convert the MLT file to the espa xml format I get this error:

Error: read_lpgs_mtl : Unsupported resampling option: CC

It turns out that the resampling option found in these metadata files is:

RESAMPLING_OPTION = "CC"

While the source code is looking for the string "CUBIC_CONVOLUTION":

if (!strcmp (tokenptr, "CUBIC_CONVOLUTION"))
    tmp_bmeta.resample_method = ESPA_CC;

When did the metadata standard change? How can I make both compatible?


texts: Shape map with a merged vector layer not working in power BI

texts:
I merged two shapefiles—one of school districts and one for intermediate school districts. I went through every imaginable step to make sure the attribute fields were the same for both layers. I then converted it into a topojson file and imported it into Power BI. I made sure the my data had matching map keys in the topojson file. When I go to create my map, K-12's show up fine, ISDs won't show up at all, even if I try to map them on their own. Is there any way to make this work?
Steps taken:

Downloaded two shape files from the Michigan Open GIS Data portal (http://gis-michigan.opendata.arcgis.com/datasets/school-districts-v17a; http://gis-michigan.opendata.arcgis.com/datasets/intermediate-school-districts-v17a/data)

Loaded them into QGIS and merged the two layers into one layer, exported as a new shapefile.

Loaded merged shapefile into mapshaper, exported as topojson file.

Chose the shape map visual in power bi, added my merged topojson file as the map. Location name is the same in the attributes table and my main data table that I'm trying to visualize, so is the district code.

When I add the school district name for K-12 districts everything shows up fine. If I try to make a map of ISD districts only, nothing shows up. It's like it doesn't recognize that they're joined even though there are linked attributes in both sources.


texts:
I am sending my modified features to database with transactional WFS. The problem is when I render my vector features on map their property name which describes geometry is changed to 'geometry' and geometry column name in database is wkb_geometry(in this case).
So every time I am sending request via WFS-T I have to set geom property name to 'wkb_geometry'

feature.set('wkb_geometry', feature.getGeometry());

, if I leave 'geometry' response will return :

<ows:ExceptionText>No such property: geometry</ows:ExceptionText>

Is there some workaround for this because I can't hardcode it this way, some times maybe geometry column name will be something else except wkb_geometry.


texts: How to reuse Thematic layer definitions in Mapinfo?

texts:
I have a map from table. Let the table be TABLE1. Now That I have created Thematic Map from Column1 of this table. Again I executed the following query

Select * from TABLE1 where TABLE1.column1 ="someVal" into TABLE2
Add Map auto  Layer TABLE2

Now that new table has been created the old thematic layer does not work. Again I create new thematic layer from TABLE2 and execute following query

Select * from TABLE1 where TABLE2.column1 = "otherVal" into TABLE2.
 Add Map auto  Layer TABLE2

Thus again thematic map I created is lost.

I have created buttons to execute these query but each time I click the button the thematic layer is lost. Is there any way that I could preserve the thematic layer and apply it to TABLE2 for all select query executed.

Actually the motive was to use thematic map that was created for TABLE1. But as the points or value change in each selection and points of thematic layer remains constant and lie in same initial coordinate the thematic layer becomes unusable.

  What is the right solution for reusing thematic layer in each Query?


texts: OSM data: "Referenced nodes" vs regular nodes

texts:
I am using osm2po software to extract nodes and edges from europe-latest.osm.pbf. When run from terminal, the logs are as follows:

So, there are over 2 billion nodes, but only (or "only") around 197 million "referenced" nodes. Now my question is, what is a referenced node?


texts: Is there a way to select polygons by geometry (right angle)?

texts:
Is there a way to select polygons by geometry?

I need to find a way to detect all the polygons that have a right angle, but I can't find a tool that can do that.

I'm using QGIS 2.10.1 and ArcGIS for desktop 10.3.
I'm open to any other free software.


texts:
Using the arcgis REST api (http://resources.arcgis.com/en/help/rest/apiref/) using the addFeatures link http://example.com//arcgis/rest/services/name/name/FeatureServer/0/addFeatures how do you add a blob. We've tried the following,

[
{
   "geometry" : {"x" : -118.15, "y" : 33.80},          "attributes" : {
     "name" : "Joe Smith",
     "title" : "Test",
     "descr": "324DFIOWNF83294"
   }
 }
]

does not work (error being Unable to complete operation,   "Setting of value for descr failed.") however

[
{
   "geometry" : {"x" : -118.15, "y" : 33.80},          "attributes" : {
     "name" : "Joe Smith",
     "title" : "Test",
     "descr": null
   }
 }
]

edit Base64:

I also tried the following base64 version as Vince suggested

VGhpcyBpcyB0ZXh0IGVuY29kZWQgaW4gYmFzZTY0Lg==

being

This is text encoded in base64.

with the following,

[
{
   "geometry" : {"x" : -118.15, "y" : 33.80},          "attributes" : {
     "name" : "Joe Smith",
     "title" : "Test",
     "descr": "VGhpcyBpcyB0ZXh0IGVuY29kZWQgaW4gYmFzZTY0Lg=="
   }
 }
]

doesn't work.

What encoding does ArcGIS's REST service expect the blob to be in.

The field types are

name - esriFieldTypeString

title - esriFieldTypeString

descr - esriFieldTypeBlob

Specifically descr is,

descr ( type: esriFieldTypeBlob , alias: desc , editable: true , nullable: true )

ArcGIS Server version is 10.3.1


texts: Projection problem with ESRI JSAPI and ArcGIS map service

texts:
I was trying to obtain a gpx file with some coordinates by drawing on an Openlayers map with an ArcGIS baseMap.

When I draw the polyline and create the gpx, if I open it on Google Earth, what I see is not what I drawed before, the line is totally different from the original and not positioned where I drawed it.
I know it's a projection problem, I've tried trasforming the geometry object from Mercator to Geographic, also getting directly the geographic coordinates from the map coordinates, but nothing.
I tried to set "spatialReference" to 4362 and then to 3857, but nothing changes.

I'm going to use that .gpx on a gps device (the next week I'll go to the Svalbard islands and I need some gps tracks to go around Longyearbyen by snowmobile, there there aren't any sign of life out the town, so I must be prepared to it), when I'll be there I'll adjust the output right for the device they will rent to me, but now I need to save on the .gpx file almost the right coordinates.

I'm getting from the map those coordinates:
lat: 61.22582068741976
lon: 4.684820015391338
when I'm expecting instead something around 78. lat and 15. lon.

This is some of the code I use to create the map (I'm not pasting the code I know it's not responsible of my problems):

    var initialExtent = new esri.geometry.Extent({"xmin":505615.5801124362,"ymin":8678955.717187276,"xmax":525935.6207525175,"ymax":8689168.654279819,"spatialReference":{"wkid":32633,"latestWkid":32633}});

    map = map = new esri.Map("map", {extent: initialExtent, logo : false});
    basemapURL = "http://geodata.npolar.no/ArcGIS/rest/services/inspire1/NP_TopoSvalbard_U33_CHL/MapServer/";
    map.addLayer(new esri.layers.ArcGISTiledMapServiceLayer(basemapURL));

Here I'm using wkid 32633 that is the default for that map, tried to change with known ones, but nothing happened.

And now the code I use to get the coordinates:

        dojo.connect(tb, "onDrawEnd", getCoords);
        function getCoords(geo) {
            var r=confirm("Salvare tracciato?");
            if (r==true) {
                geo = esri.geometry.webMercatorToGeographic(geo);
                for ( var path = 0; path < geo.paths.length; path ++ ) {
                    for ( var pt = 0; pt < geo.paths[path].length; pt++ ) {
                        tra += '<wpt lat="'+geo.paths[path][pt][1]+'" lon="'+geo.paths[path][pt][0]+'"></wpt>';
                    }
                }
            }
        }

"tra" is a variable that stores all the code I'll insert into the gpx file with an other function.
The "webMercatorToGeographic" function transform the map coordinates to geographic ones.


texts: Importing all layers available on server for WFS in PyQGIS?

texts:
I'm trying to access data from a WFS-server,
E.g. https://geo.barentswatch.no/geoserver/bw/ows?

I manage to import single layers by:

uri = "https://geo.barentswatch.no/geoserver/bw/ows?srsname=EPSG:4326&typename=bw:npdsurveyongoing_latest&version=1.0.0&request=vlayer=QgsVectorLayer"
vlayer = QgsVectorLayer(uri, "bw:npdsurveyongoing_latest", "WFS")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

However, I want to import all available layers and would like them to automatically to be named by their original typeName. Seems a bit unnecessary  to hardcode all names in beforehand accessing the server, but I have been struggling to figure out the codes for this.

Any way to do this in pyqgis?


texts: Earth Engine Image Collection Reduction Glitch

texts:
I am having a visual glitch when adding an image as a band to the map. The glitch consists of a vertical line which is not rendered at all. The more I zoom out of the map, the more this line increases in width, When zoomed in, this line decreases in width.

The problem is specifically when I reduce a collection to the last element while filtering its bounds.

var image = imageCollection3.filterBounds(Lixin);
image = image.reduce(ee.Reducer.last());
var NDVI = image.normalizedDifference(['B8_last', 'B4_last']).rename('NDVI');

When I try to do the same process by reducing to the sum or the mean, the problem does not happen.

var image = imageCollection3.filterBounds(Lixin);
image = image.reduce(ee.Reducer.mean());
var NDVI = image.normalizedDifference(['B8_mean', 'B4_mean']).rename('NDVI');

It is also worth noting I don't get any errors in the console.

This is the link to my code:
https://code.earthengine.google.com/f512a8c70b5f0d9444285cf190095f43

Is there a way to correct this?

How it looks when reducing to mean/sum:

How the band looks when reducing to last.


texts:
I am looking to build a 3D Web Mapping Application using data/projects initially created within GOCAD and SKUA. We currently use ArcServer 10. Initial looking around online has suggested that I could export the GOCAD and SKUA data as TINs, importing the data into ArcGlobe and then serving it out via ArcServer.

It seems as if this is way too easy and could not possibly handle the geological layers within the data. Does anyone know of an alternative method?


texts: How to get pure base topo map without any human made things?

texts:
I need detailed world map (preferably vector), just like Google Maps Terrain, but without roads, streets, airports etc. Just pure base map like an atlas with mountains, deserts, rivers, lakes etc. Preferably vector because I need a high detail level (small island, coastline etc.).
I've tried a lot of sources but they're not detailed enough or they contain a lot of roads, names, buildings... which I couldn't remove.

EDIT: And, why this url gives me a BLACK image? Shouldn't it be a terrain?
https://mt1.google.com/vt?lyrs=t&x=1326&y=3143&z=13


texts:
I'm new with using servers but I've been using ubuntu for more than a year now. I've tried to install Zoo-project by following the instructions on the website : http://zoo-project.org/docs/install/debian.html. At first I tought that it was easy but when I try to reach the testing page : http://127.0.0.1/cgi-bin/zoo_loader.cgi?ServiceProvider=&metapath=&Service=WPS&Request=GetCapabilities&Version=1.0.0 I've an error message saying :

" This site can’t be reached

127.0.0.1 refused to connect.
Try:
Checking the connection
Checking the proxy and the firewall "

I have the feeling that there is something I'm doing wrong. I don't know what... I've never started an server.

When I'm doing HTML, I start a server with python2 -m SimpleHTTPServer 8000 and I can see it in the port 8000 of 127.0.0.1 but when I try to use ZOO-project the instructions doesn't telling me how to start the server.

Do I need to install an apache server on my computer ? What is the step that I'm missing here ?

I don't know if I'm clear enough... I'm a GIS french guy after all lol.


texts: Enterprise Database Workflow with Online Mapping

texts:
We bought a trimble units for our operators to collect data out in the field using Collector. We currently have ArcGIS Server 10.3.1 with SQL Server 2012 setup on one of our servers. Our server manager has been federated with portal and all of our data has been versioned. All of our layers currently have the Z & M values turned on.

From what I've read and tested these values will not work for the online mapping use.

Do I have to recreate the database and each layer without Z and M Values?


texts:
I am currently trying to edit a layer in postgis by using qgis, but still out of luck. I don't know what I am doing wrong, and I guess it must be a very basic issue since I haven't found any question like this one.

The versions are:

postgresql 9.3.5 + postgis 2.1.7
qgis 2.12.2

These are the steps I have already followed. First I created a simple table with a geometry column

After that, I made a connection in qgis

Then, add it as a layer in canvas. At the end, when I start editing I see an error in the log window and I cannot figure out what is it about. There seems to be a problem in the cursor definition.


texts: Points dissappearing after pasting symbology from another layer in QGIS 3.10

texts:
I am planning to plot some graduated points, based on their value, while also keeping the scale the same for different parameters. Naturally, I imported the CSV file that included the points. Upon plotting as categorized points, I set the value for the corresponding parameters, as you can see in the screenshots;

The points were automatically classified in 8 classes for both cases and these returned me with the following plots:

Since I want to compare these two parameters on the same scale, I copied the symbology from the Sal_Sat layer and pasted it upon SS_Model-Bias then changed the SS_Model-Bias layer's value to its original parameters which looked something like this:

Applying one layer's symbology to another, made some points disappear which was visible in the earlier classifications.

Could anyone kindly tell me what is causing the points to disappear after symbology of another layer of different style is being pasted on it? Is it the right way of making two scales the same or perhaps there are other ways of doing so?
What am I doing wrong? I am using QGIS 3.16 LTR.


texts:
I have some rasters with different morphological features each (slope, exposure etc), There is a way to select points (from a POINT shapefile) by a raster characteristics? (i.e. points in area with slope >25%).

There is, moreover, the possibility to select points combining 2 or more geomorphological feature (i.e. Point in area with slope >25 & Exposure from 45° to 90°)

I am using QGIS v2.18


texts: Split a network according to a vector of nodes

texts: I would like to split a network of roads at the vertices just created in another layer. So, I have got a roads layer and a nodes layer and I would like to exploit the nodes layer in order to obtain a graph from the roads made by edges defined by those nodes.


texts: WFS geotools Hashmap can not be parsed to SimpleFeature

texts:
I'm trying to implement a WFS client in Java using the Geotools. But everytime I tried to read the actual features I get following Error

java.lang.ClassCastException: java.util.HashMap cannot be cast to org.opengis.feature.simple.SimpleFeature

It seems it is nested somewhere in

  (PullParserFeatureReader.java:84)
  Requests to other WFS services work. So it seems to be something specific to this service, which I can't figure out.

My code looks like this

    Map connectionParameters = new HashMap();
    connectionParameters.put("WFSDataStoreFactory:GET_CAPABILITIES_URL", "https://geodaten.metropoleruhr.de/poi/poi?" );
    connectionParameters.put("WSFDataStoreFactory:WFS_STRATEGY", "mapserver");
    connectionParameters.put("WFSDataStoreFactory:MAXFEATURES", "1");

    WFSDataStoreFactory dsf = new WFSDataStoreFactory();
    WFSDataStore data = dsf.createDataStore( connectionParameters );

    String typeNames[] = data.getTypeNames();
    String typeName = typeNames[0];

    try {
        SimpleFeatureSource source = data.getFeatureSource(typeName);
        SimpleFeatureCollection fc = source.getFeatures();

        while(fc.features().hasNext()){
            SimpleFeature sf = fc.features().next();
            System.out.println(sf.getAttribute("name"));
        }
    } catch (IOException ex) {
        ex.printStackTrace();
    }


texts: Geocode incredibly quickly but inaccurately

texts:
I have a Rails 4 app that needs to be able to geocode a large number of places very quickly, but not with any high degree of accuracy. Accurate to a few miles would be fine, I really just need the lat/lon of the nearest city.

Some solutions I have thought of would be to download my own database and use some sort of levenshtein distance function to pull the correct record by city/country, but ideally there would be some kind of API that geocodes innacurately, but quickly and cheaply.

I need to be able to run possibly hundreds of thousands-millions of queries a day, with sub 2ms type speed.

A "place" in my data set is user-inputed. Most likely a city/country or city/state. e.g. "Santa Barbara, CA" or "Ontario, Canada"


texts:
I am having issues with my layers. I am using QGIS v3.2, and have been stuck on this issue for a while now. I using OGA License blocks to plot certain points at sea. These license blocks follow a standard square shape e.g 4km x 4km. However, some of these blocks have multiple irregular shapes within. the Figure below identifies examples of standard squares and squares with irregular shapes within.

I need these irregular shapes deleted. The process I have used is to --> download shapefile with data. when I load the shapefile into QGIS and open the attribute table none of the column headings relate to causing the irregular shapes. I have attempted to delete the irregularities by --> toggle editing --> selection by feature and selecting the features (figure 2)

When I delete them however, although they delete and erase the irregular blocks they also erase the joining line between each 4km by 4km blocks which I need, as seen in figure 3

Does anyone know how to solve the issue of deleting the irregular shapes within the blocks whilst at the same time not deleting the boarder to each 4km by 4km block?


texts: Counting overlaps on large line dataset using FME runs into memory issues?

texts:
I have a large dataset of lines (ca 400 000). Each line represents a calculated route on the same road network.

My goal is to summarise a total number of calculated trips that pass each part of the road network, ie. a line on line overlay resulting in an overlap count on each segment.

I've tried running it in the line on line overlayer in FME, but I keep running into memory issues. I've also tried intersecting in Spatialite, but have memory issues there as well.
I'm stuck on a 32 bit computer at work, so if the answer is "that dataset should be fine on 64 bit with enough RAM" that could be a viable answer.
But something also tells me there should be a more elegant way to go about it...

So, in short, what's an efficient way of counting overlaps/intersecting in a large line dataset? I have access to and mainly work with FME, QGIS and R.


texts:
I am using pgAdmin (v.1.22) for my postgres(9.5). I am still fairly new to postgres. I am trying to perform a union(intersection) of multiple layers where there is both complete and no overlap between layers. I need all linework and attributes to come through to the final product. I am playing around with st_union and st_intersection but can't quite get to the final layer I need.

DROP TABLE IF EXISTS TEST;
CREATE TABLE TEST AS
SELECT
 (st_dump(st_union(a.shape, b.shape))).geom geom, a.status, b.code
FROM
 p760_passive_lb a, p760_fmu b
WHERE
 st_intersects(a.shape, b.shape)
;
 ALTER TABLE TEST ADD COLUMN id SERIAL;
ALTER TABLE TEST ADD PRIMARY KEY (id);
 CREATE INDEX idx_TEST_geom ON TEST USING gist(geom);

Where am I going wrong?

This is similar to:
Is there a union function for multiple layers comparable to ArcGIS in open source GIS like QGIS?

So, to help illustrate, I have attached an image of what is going on. The top two layers are my inputs. The bottom left is the output of the above code. Two things of note: (1) the features in the western portion are not visible, and (2) when I use the identify tool (bottom right) in QGIS on any ONE area in the result, the whole area gets chosen and multiple features are selected. Strangely, the linework in the western section is viewable with the identify tool. This leads me to think that the layers are just stacked on top of each other and not forming one geometry.


texts: Updating field to give count of points in polygon using STIntersects?

texts:
I have a points layer (dbo.ptLayer)

Around 1M points
Spatial Geometry Type (dbo.ptLayer.geom)
No spatial index just yet, but will create one once data gathering complete.

I have a polygon layer (dbo.polygonLayer)

Around 500 polygons.
Spatial Geometry Type (dbo.polygonLayer.geom)

Both have fields called ID.

How do I populate an empty integer field in the polygon layer, with a count of the total number of points within each polygon?

Although I have access to other software products, I am interested to learn what can be done purely within SQL and SQL Server.

I believe I should be making use of STIntersects but would like to know what is the best way of doing an update to populate this field.


texts: Speed up plug-in datasource in ArcObjects

texts:
Are there any tricks involved in getting plug-in data sources to run at higher speeds, similar to the built in types?  I've noticed that a simple plugin-data source will run many times slower than the equivelent data loaded from a shapefile or personal geodatabase, and the extra time seems to be mainly in the rendering phase.

Are there any tricks available to improve this in ArcObjects?


texts:
I am trying to update a field with a new value as I updated a Domain for a web purpose.  I am new to this.

I want be able to choose the layer and then the field to update

M=Municipal
S=State
C=County
P=Private

Here is the script so far.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Ownership Domain Update.py
# Created on: 2016-11-15 11:57:22.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: Ownership Domain Update <FClass>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Script arguments
FClass = arcpy.GetParameterAsText(0)
FieldName = arcpy.GetParameterAsText(1)

# Process: Select Municipal
arcpy.SelectLayerByAttribute_management(FClass, "NEW_SELECTION", FieldName  = 'M')

# Process: Calculate Municipal
arcpy.CalculateFieldName_management(FClass, FieldName , "\"Municipal\"", "VB", "")

# Process: Select State
arcpy.SelectLayerByAttribute_management(FClass, "NEW_SELECTION", FieldName  = 'S')

# Process: Calculate State
arcpy.CalculateFieldName_management(FClass, FieldName , "\"State\"", "VB", "")

# Process: Select County
arcpy.SelectLayerByAttribute_management(FClass, "NEW_SELECTION", FieldName  = 'C')
# Process: Calculate County
arcpy.CalculateFieldName_management(FClass, FieldName , "\"County\"", "VB", "")

# Process: Select Private
arcpy.SelectLayerByAttribute_management(FClass, "NEW_SELECTION", FieldName  = 'P')

# Process: Calculate Private
arcpy.CalculateFieldName_management(FClass, FieldName , "\"Private\"", "VB", "")

# Process: Clear
arcpy.SelectLayerByAttribute_management(FClass, "CLEAR_SELECTION", FieldName  = 'P')

And the Error:

ERROR:
raceback (most recent call last):
  File "G:\_LSSE GIS Toolbox\_LSSE Toolbox\Ownership Domain Update.py", line 18, in <module>
    arcpy.SelectLayerByAttribute_management(FClass, "NEW_SELECTION", FieldName  = 'M')
TypeError: SelectLayerByAttribute() got an unexpected keyword argument 'FieldName'

Failed to execute (Script).
Failed at Tue Nov 15 14:14:33 2016 (Elapsed Time: 0.06 seconds)


texts: Use arcpy to query by date field in ArcGIS for Desktop?

texts:
I have an ArcGIS table in a file geodatabase that contains a date field. I would like to write a python (arcpy) function that creates and returns a table view for a date-range given a starting and/or ending date.

Here's what I have so far:

def get_date_range_table_view(grid_table, start_date=datetime.datetime(1, 1, 1),
                              end_date=datetime.datetime(9999, 12, 31)):

    min_date = '"DateTime"' + " >= date '" + str(start_date) + "'"
    max_date = '"DateTime"' + " <= date '" + str(end_date) + "'"
    expression = " AND ".join([min_date, max_date])

    return arcpy.MakeTableView_management(grid_table, "date_range_table_view",
                                          expression)

The expression that uses 'start_date' works great and returns a table view that can be used with a SearchCursor. The trouble that I'm having is getting the 'end_date' logic to return a similar table view. Although the logic is identical, except for the '>=' symbol, the table view that is returned has no readable rows.

I found a similar question here:
Arcpy cursors, WHERE clauses, and date/time fields

Which directs the readers to this resource (which has information about querying Date fields in ArcGIS 10.2 based on database type):
SQL reference for query expressions used in ArcGIS


texts: Sum the field of a relationship table into the source table

texts:
I have created a relationship between two feature classes(one to many) and i want to have the sum of a field in tableB(many) into a cell in TableA (one), TableA is states and tableB is roads within the state, so i want the sum of roads in each state(in a field in tableA).
I have created one to many relationship class between two feature classes (not shapefile) within the same database, these feature classes have attribute tables (A, one and B, many). So i want sum of a field(colunm) from tableB into tableA.
I hope it is clear know.


texts: Join Multiple Tables to a Spatial Map without using Database

texts:
I work with Spatial dataset in QGIS and as usual, it is linked to an attribute table.

What is now unusual to me is that the attributes data is now growing and my client doesn't want the data merged on one table. He needs different tables as the data comes in.

So, I want to know how I can have the spatial data linked to multiple tables without necessarily using a database engine. Just the normal attributes tables in QGIS linked to one spatial dataset?

I need a way to do this without using database (such SQLite, PostGIS, MySQL Spatial etc)?

If QGIS don't support this by default, any suggestion on any Open Source GIS plugin/software is also well come, but no database.

How do I link more than one Attribute tables to Spatial data on my map in QGIS (no database)?


texts: I am designing my map with tilemill. Now I am struggling with data visualization. I would like to visualize the road name like the 1st image instead of the second image. I have loaded the data source of tilemill from shapefile which contain Road name that field I want to visualize. Can anybody help me on this issue, whether it is possible or not to make the annotation from data?


texts: Reverse geocoding to get city, state, country data separately?

texts:
I have longitude and latitude data. I need to reverse geocode to get city, state, country data separately.

I have checked google service.

Google provides complete address to a micro level. Separation into city, state, country seems not possible.

Any suggestions?


texts: Failure connecting to API - Geocoding using MMQGIS

texts:
I'm very new to GIS and am trying to geocode a list of addresses of some OpenData in QGIS using MMQGIS....however when I go to run MMQGIS Geocode function with my .csv. I keep getting this error and similar ones.

Failure connecting to API: <urlopen error [SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1056)>

I've read the different forums and updated my Shell Profile and Install Certificates in Python. If its relevant I use MacOS and QGIS 3.10 and have Python 3.8 Installed..


texts:
Could I load spatial data to Oracle when table structure was created before this step? Specifically I have problem with default name of attribute OGR_FID . If my primary key has different name, I can't load data.

ogr2ogr -append -f "OCI" OCI:abcd/abcd@(DESCRIPTION=(ADDRESS=(PROTOCOL=tcp)(HOST=a000000)(PORT=1521))(CONNECT_DATA=(SID=hop))) PG:"host=w111111 user=post dbname=gis password=gres" -sql "select row_number() OVER (ORDER BY null) fid,typppd,geom from zab.table" -lco MULTI_LOAD=yes -lco GEOMETRY_NAME=sdo_geom -nln HRANICE

I found Layer Creation Options -lco FID=id in PostGIS driver implementation, but OCI ignore this in my opinion.


texts: Where can I obtain Water Virtual Flow Segments for Minnesota?

texts:
This is a cross-border task.  I already have this data for Ontario, Canada but I also need to complete the model for a large lake that is essentially split in half between the two countries.

Definition of Water Virtual Flow Segments:
Depicts the flow of water from one stream to another through a lake or wetland.

Here is an example of the data I have (Ontario Ministry of Natural Resources):
I am looking to find a comparable Minnesota dataset.


texts:
I need to export the current OpenLayers view into a canvas to finally save it as a image.

My problem is, that I need an image of the current user view with some interactive elements.

I saw this working example:
http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/examples/exportMapCanvas.html

There is a little addon that is not OpenLayers standard, but add a little export functionality:
http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/examples/exportMapCanvas.html

but this does no work on current open layers (2.13) any more and I don't get which version that is.

The Constant OpenLayers.Layer.Grid.ONECANVASPERLAYER does not exist (any more ?) And the option "useCanvas" is deprecated as well.

This is the broken code (summarized):

        if (layer instanceof OpenLayers.Layer.Grid){
            this.drawLayer(layer.canvas);
        }

A Grid layer has no attribute "canvas" as far as I can see it.


texts: How do you delete Survey123 attachments from a survey?

texts:
There seems to be no documentation on how to delete Survey123 attachments. Is it actually feasible?

Is deleting the complete feature layer record required/necessary? (and if so, how can one be sure that the attachment was deleted, since the link to it was deleted?!)

Can the survey record be retained, but the attachment & relationship deleted?

Edit: My surveys and their attachements are stored in an ArcGIS Online Feature Service.


texts: About a widget to get coordinates in Web API?

texts:
I need to know if there is a widget that lets me draw a geometric figure in the map, and get the coordinates of the selected area.
Any similar widget is welcome.


texts: How to get information about layers using Arcpy

texts:
I have this layer structure in ArcMap:

"Localities, Gazeteers, Bases, Regions, and Territorial Claims are all Group Layers. All others under these are layers. I have a script that uses the Arcpy Describe method to go under each layer and get its full directory path. Here is the script:

import arcpy
mxd=arcpy.mapping.MapDocument("CURRENT")
dataframes=arcpy.mapping.ListDataFrames(mxd)
for dataframe in dataframes:
    for lyr in arcpy.mapping.ListLayers(mxd, "", dataframe):
        if not lyr.isGroupLayer:
            print(arcpy.Describe(lyr).catalogPath)

This would return something like: C:/Folder/Subfolder/layer.shp
Now, I am trying to figure out a way to get the path of each layer relative to the map. For example, for layer "features_multipoint", I need the script to report something like this:

MapDocument/Localities/Gazeteers/features_multipoint

The Describe method seems to not provide such capabilities.


texts: Understanding curvature filter of QGIS raster terrain analysis?

texts:
I have read the source code of several QGis-1.7.4 raster filters computing slope, aspect and curvature.

There is a formula in the filter computing the total curvature that puzzles me.

The source file is in the current version of QGis, with the following path:

qgis-1.7.4/src/analysis/raster/qgstotalcurvaturefilter.cpp

The purpose of this filter is to compute the total curvature of the surface in a nine-cells window. The function code is as follow:

float QgsTotalCurvatureFilter::processNineCellWindow(
   float* x11, float* x21, float* x31,
   float* x12, float* x22, float* x32,
   float* x13, float* x23, float* x33 ) {

  ... some code deleted ...

  double dxx = ( *x32 - 2 * *x22 + *x12 ) / ( 1 );
  double dyy = ( -*x11 + *x31 + *x13 - *x33 ) / ( 4 * cellSizeAvg*cellSizeAvg );
  double dxy = ( *x21 - 2 * *x22 + *x23 ) / ( 1 );

  return dxx*dxx + 2*dxy*dxy + dyy*dyy;
}

I'm ok with the "dxx" formula, and with the return expression. But I think that the "dyy" and the "dxy" formulas are inverted: this makes the total result asymmetric regarding x and y dimensions.

Am I missing something or should I replace the double derivatives expressions by:

  double dxx = ( *x32 - 2 * *x22 + *x12 ) / ( 1 ); // unchanged
  // inversion of the two following:
  double dxy = ( -*x11 + *x31 + *x13 - *x33 ) / ( 4 * cellSizeAvg*cellSizeAvg );
  double dyy = ( *x21 - 2 * *x22 + *x23 ) / ( 1 );
  return dxx*dxx + 2*dxy*dxy + dyy*dyy; // unchanged

Could you tell me your opinion on these formulas, if they are incorrect as I thought or if I am wrong? In this last case, do you know why the formulas have to be asymmetric regarding x and y?


texts: how to select line segments not overlapped by other layer's lines

texts:
I have two shapefiles of forest roads (mostly dirt) from different sources, covering the same area.  The goal is to build one shapefile of roads that is the composite of the other two.

Layer A is less complete, and Layer B is more complete, i.e. if a dirt road on A is a mile long, the same road may exist on B but is a mile and a half long, and the first mile mostly overlaps with the A line but not exactly.

The trivial answer here is "just take the B road and ignore the A road" but there are cases that are reversed from the example, i.e. there are some cases where A has more detail / more roads than B.

how would you do overlap detection?  Buffer each layer and look for overlaps?
how would you select only the portions of lines that are not in the overlapped region, for example, the last half mile of the B road from the example?


texts: How to iterate through the attributes of a feature in Leaflet?

texts:
I receive polygons via the ESRI REST API with Leaflet. I want to iterate through all attributes of a feature in order to create a custom pop-up for it.

My code looks like this:

myLayer.on('click', function (evt)
{
   feature = evt.layer.feature;
   //Here I want to iterate ...
}

I know that I can access the attributes like this:

feature.properties.myAttribute

But how can I iterate through all attributes?


texts: How to calculate the average point spacing for a layer of points in QGIS?

texts:
I have a set of points and want to calculate the average point spacing for further analysis.

ArcMap has the 3D Analyst Point File Information for this but I'm searching for an alternative in QGIS. Any ideas?


texts: Changing unit of measure from degrees to meters in QGIS?

texts:
I am running QGIS 1.7.0 on windows, and I donwloaded a shapefile of census block group boundaries from TIGER.  This file appears to be in the WGS 84 coordinate system, and when I open it in QGIS the unit of measure is degrees.

I tried re-projecting the file into the NAD 83 coordinate system, and saving it as a new shapefile, but when I open the new file it still displays a ruler in degrees.  If I try to make a buffer, it interprets buffer distance as degrees.

Assuming I'm starting with a shapefile from TIGER, what's the best way to convert my unit of measure from degrees to meters?

For posterity, it seems like the following workflow make sense when dealing with US census TIGER shapefiles:

Download the relevant shapefile
Open it in QGIS
Click the Settings menu and select project properties
Select a projected coordinate system, like NAD83 / UTM zone 19N
Click the box Enable 'on the fly' CRS transformation
Click ok
Click the 'zoom full' button (one of the magnifying glasses)


texts: CartoDB map not viewable by anyone, despite being public

texts:
I have a map created with CartoDB, but when I give the bit.ly link to others, or view it in incognito, I see the message:

"Looks like this map no longer exists or is set as private."

This is the link to the map: http://bit.ly/1OwBmCB

The visibility is set to public.

How can I make it so that others can see this map?


texts: Permanently changing color of tool/output in ArcGIS ModelBuilder?

texts:
I have a pretty complex model in ModelBuilder, and I'd like to custom-color certain tools and outputs. I'm not so fond of the orange-green-blue scheme that is there by default.

However, when I change the color of one of the blocks using 'display properties,' the color is changed back to the default whenever I validate the model.

Is there any way to make the color change permanent?


texts: ArcGIS importing from file geodatabase A to B without stopping service

texts:
I have file geodatabase B which is use by arcgis server as a data store

and on a schedule interval I need to import new data from file geodatabase A.

geodatabase need to be wiped and accept the data from geodatabase A.

May i know if there is a arcgis tool for this?

This process must not stop the arcgis server/services


texts: Need help to convert 9.3 code for registering .dll to 10.1

texts:
This is my code to register my DLL (I'm using ArcMap 9.3)

Private Sub regDLL(ByVal dllName As String)
   Dim asm As System.Reflection.Assembly = System.Reflection.Assembly.LoadFile(dllName)
   Dim regAsm As New RegistrationServices()
   Dim bResult As Boolean = regAsm.RegisterAssembly(asm, AssemblyRegistrationFlags.SetCodeBase)
End Sub

This same code won't work in ArcGIS 10.1 since there is esriRegAsm.

Does any one know how to convert it so that it will work?

I'm using Visual Studio 10.1.


texts: r.reclass slope interpolation error on mac

texts:
Using this video to create a slope analysis for a land I am designing the following error message occurred after I tried to run the r.reclass add-on on a before interpolated file.
These are the parameters I used:
0 thru 4.398705 = 1
4.398705 thru 4.763641 = 2
4.763641 thru 6.340191 = 3
6.340191 thru 9.462322 = 4
9.462322 thru 80 = 5

Could there be something wrong with the way I wrote the parameters? And can someone please help me with the error message below? Your help is much obliged.
-----------------BEGINNING ERROR --------------
Traceback (most recent call last):
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 545, in writeCommand
output.write(command.encode('utf8') + '\n')
TypeError: can't concat str to bytes

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 434, in processAlgorithm
Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 375, in executeGrass
command, grassenv = Grass7Utils.prepareGrassExecution(commands)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 362, in prepareGrassExecution
Grass7Utils.createGrassBatchJobFileFromGrassCommands(commands)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 278, in createGrassBatchJobFileFromGrassCommands
Grass7Utils.writeCommand(fout, command)
File "/Applications/QGIS-LTR.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 548, in writeCommand
output.write(command + '\n')
UnicodeEncodeError: 'ascii' codec can't encode character '\xd6' in position 151: ordinal not in range(128)

Execution failed after 0.53 seconds
--------------END---------------


texts: What is the practical limit for lat/long measurement accuracy?

texts:
There are several systems for finding lat/long coordinates. This includes GPS satellites, DoD TRANSIT, Russian GLONASS, Galileo, differential GPS, LORAN, land-based systems, and more. Maybe even using laser interferometry directly from Greenwich.

I want to know, using all of this technology, what is the actual maximum accuracy which a location on Earth is measurable?

This is not a question about technology or measurement. The theoretical limit is nanometers because of lasers. But the actual limit will be worse because of the physical characteristics of Earth, it is constantly moving.

So what is is actual best measurement we can make?


texts:
I am a self-taught QGIS user with low computer skills.

I run QGIS 2.10 on Ubuntu 15.10. I am trying to use Alexander Hajmal's tile_writer script (http://alephnull.net/software/gis/tile_writer.shtml). However, when I attempt to execute the script I get this message:

Traceback (most recent call last):
  File "", line 1, in
  File "/home/brian/tile_writer.0.2.1/tile_writer.py", line 64, in
    from globalmercator import GlobalMercator
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 478, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named globalmercator

I have globalmercator.py on my home directory along with tilewriter.py, in a folder named "tile_writer.0.2.1". When the script didn't work I tried to add the following two lines to my etc/environment file:

PYTHONPATH="${PYTHONPATH}:/home/brian/tile_writer.0.2.1/globalmercator.py"
export PYTHONPATH

Unfortunately the results were still the same.

How do I properly set the path, or whatever else I am supposed to do?


texts:
I'm using OpenLayers 5 to show a map and add markers to it.
I can see the map without any problems but only the last marker is shown. I saw an answer about using .clone() method on the marker that's being added to solve this issue(apparently markers were being recycled) bu i do not know how to do this in OpenLayers 5. Here is my code:

var map = new ol.Map({
        target: 'map',
        layers: [
            new ol.layer.Tile({
                source: new ol.source.OSM()
            })
        ],
        view: new ol.View({
            projection:"EPSG:4326",
            center: [24.038085,31.780256],
            zoom: 6,
            minzoom: 6,
            maxzoom: 12
        })
    });

I add markers like this:

for (var i = 0; i < 10; i++) {
                    var marker = new ol.Overlay({
                        position: [i, i],
                        positioning: 'center-center',
                        element: document.getElementById('marker'),
                        stopEvent: false
                    });
                    map.addOverlay(marker);

my marker element is defined as:


texts: Replace Polyline Geometry on Zero-Length Polyline

texts:
I have an app that occasionally generates zero length polylines with identical start and end vertices.  I need to modify the geometry for these features so that the end vertex is slightly different than the origin vertex.  I am using an arcpy.da.UpdateCursor to replace the geometry of the features, but am getting an AttributeError: 'Polyline' object has no attribute 'setValue' on the second to last line.  Any thoughts on how to resolve?

import arcpy

oilZoneFC = "T:\\Temp\\processing.gdb\\TEST"

cursor = arcpy.da.UpdateCursor(oilZoneFC, ["SHAPE@"])
for row in cursor:

    length = row[0].getLength('PLANAR', 'METERS')
    arcpy.AddMessage("Length = " + str(length) + " meters")
    if length == 0.0:
       arcpy.AddMessage("ZERO LENGTH FEATURE TO BE MODIFIED")
       for part in row[0]:
           origX = part[0].X
           origY = part[0].Y
           newX = origX + 0.0001
           newY = origY + 0.0001

           newCoord = arcpy.Array([arcpy.Point(origX, origY), arcpy.Point(newX, newY)])

           newLine = arcpy.Polyline(newCoord)
           row[0].setValue("SHAPE", newLine)
           cursor.updateRow(row)


texts: Error message while trying to load spreadsheet: Failed to connect to database. Class not registered

texts: I am trying to load an Excel spreadsheet into ArcMap 10.6 in Windows 10. But it will not load and this is the error message:


texts: Processing LANDSAT 8 data in R (using packages Raster, RSToolbox, etc.) in order to calculate NDVI and land cover classification?

texts:
I am analyzing several LANDSAT scenes from the LANDSAT 8 OLI mission (18 in total, all from Jan 2017 - cloud cover < 20%).

I now want to calculate the NDVI and do some land cover classification on these scenes using R and its RS packages, but first need to apply some corrections to these scenes.

If I'm correct, I need to (with examples of codes provided from RStooolbox and raster packages):

1. Do the DN calibration

p198r053_rad<-radCor(p198r053, metaData = metap198r053, method = "rad")

2. Apply an atmospheric correction

p198r053_haze<-estimateHaze(p198r053, darkProp = 0.01, hazeBand= c("B1_dn", "B2_dn", "B3_dn", "B4_dn"))
p198r053_sdos<-radCor(p198r053, metaData = metap198r053, hazeValues= p198r053_haze, hazeBand= c("B1_dn", "B2_dn", "B3_dn", "B4_dn"), method="sdos")

3. Apply a topographic correction (C-correction)

using a DEM and RStoolbox::topCor

4. Mask the clouds using the QA layer

using raster::mask()

5. Check the geolocation of each scene

I don't know how to do so yet (see questions below).

6. Calculate the NDVI for each scene

using

ndvi_formula<-function(nir,red) {
  (nir - red)/(nir + red)
}

(i.e. bands 5 and 4 of L8 data - but from which layer?? (see question below))

7. Classify images according to land cover type

(perhaps this should be done with the Semi-Automatic Classification tool in QGIS?)

8. Check if all scenes have the same projection and reproject if necessary

Using raster::projection() and raster::projectRaster()

9. Create a mosaic of all scenes for NDVI and land cover layers

Using raster::merge()

10. Crop mosaïcs to match the country's borders

Using a shapefile of the country and raster::mask()

(I have based my approach on the one given by Wegmann et al. (2016) in "Remote Sensing and GIS for Ecologists using Open Source Software")

Are these steps correct and in the right order?

Also, more specifically:

Steps 1-4 seem to create different rasters for each scene, how do I merge these rasters together to have a scene with all corrections from which to calculate NDVI and do the classification?
How do I correct the geolocation of my scenes?


texts: Distributed WMS (and TMS, optionally) server architecture, GPU rendering

texts:
What if we could use the power of dozens of servers to process a lot of incoming WMS requests? And more, using the power of GPU's to do the job of styling the maps for scalability and speed!

The theme of my senior research project is exactly this, distributing the WMS service (and probably, TMS) to reach a maximum horizontal processing of GIS data.
I'm very excited about this, and I would like to hear from you, the professionals in this field, ideas and criticisms regarding this approach.


texts: I am geocoding addresses using the ArcGIS Online Geocoder and Street and Zip Code. In some cases, the geocoder has read the zip code as 4 digits and in other as 5. I got a 99% hit rate, but the addresses not are geocoded correctly as a result of this error. How do I fix this? Maybe add a field for city?


texts:
I am trying to Georeference CAD Drawing in ArcMap  by adding the Control Points. But it is taking a long time to load the CAD file. I could not do snapping and every time I hover over the mouse, it will load try to load. I tried to create spatial index (under properties) but it does not help much. I realized that it is slower when i have the base map i want the to in the georeference the cad file to. I only have one CAD Polylines layer turned on.

It is hard for me to Georefernce the CAD file in this way. Is there any way to make ArcMap read/load the CAD file faster?


texts: Output raster extents when using the Spatial Analyst Extract By Mask tool to clip a raster

texts:
I am using ArcGIS 9.2 and I have a simple python script that does the following :

Loop over a polygon feature class (shown in red in the image
below)
For each record, create a feature layer to isolate the
    polygon for the row
Clip a raster (floating point DEM) with the
        geometry of the polygon feature layer using the Extract By Mask tool

The script works great; however, the clipped raster extents match those of the bounding rectangle of the ENTIRE polygon feature class (all rows). The pixels outside of the individual polygon used as the in_mask_data parameter are NoData values (shown in blue in the image below), so they can be easily dealt with.

I would like to know if it is possible to limit the output raster extents to those of the bounding rectangle of the individual polygon used as the in_mask_data parameter? I looked into the Output Extent environment setting and the MINOF setting is the default and it is supposed to be the intersection of inputs - the input raster and the feature layer in my case).


texts: Why is ModelBuilder not outputting Make Feature Layer results?

texts: This model is giving me some trouble. When I run it the is not input  into the map and when I do add the data in it is not clipped to the boundary. Any ideas? I have add model to map check in geoprocessing options.


texts: How to download Google location timeline raw data (red dots)?

texts: So I downloaded Google location timeline history. However, this data is represented as a polyline with just about a dozen time marks (Figure 2). At the same time when I look at the Google location timeline history in Google Maps, I can see so-called "raw data" represented with red dots, hundreds of them (Figure 1). So, the question is: how to download this Google timeline history point data or at least extract more time marks from a polyline of my path I downloaded?


texts: How to slice an image in gee into multiple images

texts: I want to slice a big image into images with shape 256x256xchannels, one of my channels show the class of a pixel, which can be 0, 1 or 2. My goal with this slice is to throw out the images that only have pixels with the class 0, so that when I export those images I would spend less space in memory. How can I do that in the google earth engine?


texts: Calculating Canopy Cover percentage in buffer using ArcMap?

texts:
I am trying to calculate percent canopy cover within a buffer for each site I have. Here are the steps I took.

I used the tabulate area function with the canopy cover layer to get the area for each location
I added up rows 3-74 to get total canopy (I excluded row 2 from this part since the value is zero which would mean that there isn't any canopy from what I interpreted). To get the total area I added row 2 to the total canopy number.
I divided total canopy by the total area to get the percent canopy for each site.

By using the tabulate area function I think the numbers for each row should be meters squared.

Is this correct?


texts: Using "Project (Data management)" in ModelBuilder using ArcMap

texts:
I have tried to change the coordinate system of a shapefile with point data with Project in ArcMap 10.8's ModelBuilder but nothing seems to work. I keep getting Error 99999.
I have tried the following;
Solution 1: Shape -> project
Solution 2: Shape -> feature to feature -> project
Alternative solution: CSV -> Make XY event layer -> Feature class to feature class -> project
When I try to do it manually on the same files it works just fine


texts: Creating rectangle on centroid point by map click event ol5

texts:
I am trying to develop a web portal using ol5 API where user click on map and then it takes map click event coordinates as centroid of rectangle and create a rectangle of particular height and width.

Here is the code to create square

image: new ol.style.RegularShape({
        fill: fill,
        points: 4,
        radius: 30,
        angle: Math.PI / 4
    })

But I want to create rectangle with center point, height and width.
How to achieve this?


texts: Create near table to show the distance to every feature chosen in arcgis toolbox

texts:
I'm scripting a toolbox in which the user chooses a feature containing multiple points to create a distance table from other selected features in the list.

The following script for every loop overwrites the distance. I want every loop output to be written in a new column for every feature.

Any idea?

import arcpy

Stations = arcpy.GetParameterAsText(0)
Dest = arcpy.GetParameterAsText(1)
Out = arcpy.GetParameterAsText(2)

in_features = Stations
near_features = Dest
out_table = Out

for near_features in Dest.split(';'):

out=arcpy.GenerateNearTable_analysis(in_features, near_features, out_table)


texts: Adding custom projection to epsg.properties?

texts:
I have created an entry at www.spatialreference.org of a customized .prj file.

http://www.spatialreference.org/ref/sr-org/7973/

I'm trying to add the generated geoserver entry to the EPSG.Properties file.  Upon a restart, Geoserver does not recognize my custom projection that I have attempted to add.  Reviewing the logs is not indicating any issues, its as if I didn't change the EPSG.properties file at all.  I have confirmed, in the log, that I'm using the right EPSG.properties file.

The mystery has deepened, but I think I resolved my immediate problem through trial and error.  My EPSG.properties file had another entry in it that I had successfully added months ago, and it was loading fine.  So, copied its entry, and then replaced values and parameters from the 7973 projection I wanted.  Change, Restart, Check it.  Eventually, I was able to change all the parameters, and changed the description, and it shows up in my SRS List.

Where the mystery deepens, I then tried to take out the junk entries between my two good entries, and it broke.

Here are my EPSG Entries:

/This one works/
6948=PROJCS["NAD_1927_StatePlane_Wisconsin_South_FIPS_4803",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",2000040.0],PARAMETER["False_Northing",-780.0],PARAMETER["Central_Meridian",-90.0],PARAMETER["Standard_Parallel_1",42.73333333333333],PARAMETER["Standard_Parallel_2",44.06666666666667],PARAMETER["Latitude_Of_Origin",42.0],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG","6948"]]

/*This one does not work, and originated from www.spatialreference.org */
7973=PROJCS["NAD_1983_HARN_UTM_Zone_18N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",1640416.666666667],PARAMETER["central_meridian",-75.0],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG", "7973"]]

/*This was my copied 6948, which I modified one parameter at a time, it worked, but then removing 7973 above broke it! */
96948=PROJCS["NAD_1983_HARN_UTM_ZONE_18N_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-75.0],PARAMETER["Latitude_Of_Origin",0],PARAMETER["scale_factor",0.9996],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG","96948"]]

/*So, then I added it a second time, as EPSG = 97973 */
87973=PROJCS["PHI_NAD83_UTM_ZONE18N_FEET",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-75.0],PARAMETER["Latitude_Of_Origin",0],PARAMETER["scale_factor", .9996],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG","87973"]]

/*Again, deleting the 7973 and 96948 entries broke it again!  Not only that, but restoring the old file was also broke!  But rather than re-add, I renumbered 97973 to 87973 and it now works.  I'm leaving the junk entries in there! */

I don't know what to make of it, but its in there now.


texts: What options are there for geolocation-aware SQL databses?

texts:
I need to implement an SQL database that can make queries using geo-location  based logic.

e.g.
select rows within 500m of coordinate x, y

What options are available to do this with sqlite, mySQL, postgreSQL, mariaDB, or something similar?


texts: Connecting local QGIS to Docker container?

texts:
I have a Docker container from PostgreSQL image installed.

How do I connect my local QGIS app to my PostGIS database to view the map?


texts: ArcGIS 10.3 Feature To Polygon tool not working over dateline

texts: I am working in ArcGIS for Desktop 10.3. I have a set of closed polylines that are building several distinct areas in the Pacific. The goal is to make polygons from these lines. The Feature To Polygon tool in the Data Management Tools works fine for all closed lines not crossing the dateline. However, it is not working when lines are crossing the dateline (180 degree). I get empty outputs with the warning: 'error 000117'. I tried working in different projections (GCS WGS84, PlateCaree-centered at 180 degree) and I changed the working extent to my needs, as recommended in the ArcGIS Help to no avail. I am aware that I can edit new polygons and trace the lines to achieve my goal. Is there an automated way I can change a set of closed lines, that cross 180 degree, into a polygon using ArcGIS?


texts:
I've built GDAL so many times but this time I want to add Geospatial PDF reading capability. First I built PoDoFo, but found out that Poppler is better. I found out that Win-KDE contains prebuilt binaries and libraries of Poppler. I installed KDE and modified nmake.opt as follows:

POPPLER_ENABLED = YES
POPPLER_CFLAGS = -Ie:/kde/include -Ie:/kde/include/poppler
POPPLER_HAS_OPTCONTENT = YES
POPPLER_BASE_STREAM_HAS_TWO_ARGS = YES
POPPLER_0_20_OR_LATER = YES
POPPLER_0_23_OR_LATER = YES
POPPLER_LIBS = D:\KDE\lib\poppler-cpp.lib D:\KDE\lib\freetype.lib D:\KDE\lib\lcms-1.lib advapi32.lib gdi32.lib

nmake (vc120 - VS2013) resulted without any errors. Then I copied following files to my bin dir:

freetype.dll
iconv.dll
JPEG.dll
lcms-1.dll
lcms2.dll (Dependency Walker found dependency)
openjpeg.dll
pdftoppm.exe
poppler-cpp.dll
tiff3.dll
zlib.dll

Now GDAL is unable to open PDF files. gdalinfo results in "unknown format". What do you think is wrong?


texts: Specifying distance between arrows in line SLD in GeoServer?

texts:
I have designed an SLD for a line format shapefile in GeoServer, in which I am using an arrow symbol as shown in figure.

The code is:

<?xml version="1.0" encoding="UTF-8"?>
<sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
    <sld:NamedLayer>
        <sld:Name>Blue arrows</sld:Name>
        <sld:UserStyle>
            <sld:Name>Blue arrows</sld:Name>
            <sld:FeatureTypeStyle>
                <sld:Name>name</sld:Name>
                <sld:Rule>
                    <sld:Title>A blue line with end arrows</sld:Title>
                    <sld:LineSymbolizer>
                        <sld:Stroke>
                            <sld:CssParameter name="stroke">#73ed5d</sld:CssParameter>
                            <sld:CssParameter name="stroke-width">3.26</sld:CssParameter>
                            <sld:CssParameter name="stroke-linejoin">bevel</sld:CssParameter>
                            <sld:CssParameter name="stroke-linecap">square</sld:CssParameter>
                        </sld:Stroke>
                    </sld:LineSymbolizer>
                    <sld:LineSymbolizer>
                        <!--<sld:VendorOption name="placement">centralPoint</sld:VendorOption>-->
                        <sld:Stroke>
                            <sld:GraphicStroke>
                                <sld:Graphic>
                                    <sld:Mark>
                                        <sld:WellKnownName>shape://oarrow</sld:WellKnownName>
                                        <sld:Fill>
                                            <sld:CssParameter name="fill">#0000FF</sld:CssParameter>
                                            <sld:CssParameter name="fill-opacity"
                                                >0.5</sld:CssParameter>
                                        </sld:Fill>
                                        <sld:Stroke>
                                            <sld:CssParameter name="stroke"
                                                >#0000FF</sld:CssParameter>
                                            <sld:CssParameter name="stroke-width"
                                                >1</sld:CssParameter>
                                        </sld:Stroke>
                                    </sld:Mark>
                                    <Size>20</Size>
                                </sld:Graphic>
                            </sld:GraphicStroke>
                        </sld:Stroke>
                    </sld:LineSymbolizer>
                </sld:Rule>
            </sld:FeatureTypeStyle>
        </sld:UserStyle>
    </sld:NamedLayer>
</sld:StyledLayerDescriptor>

I have to specify distance between this arrows.