texts: slmtest in R - Error non-conformable arguments

texts:
I'm running panel regression models in R with the specified model:
  tot_crime <- log_rt_crime ~ ndvi + log_disad_inde + log_density

I've successfully run:

non-spatial OLS pooled model using plm
pooled_nosp_crime <- plm(tot_crime, data = final_long_panel, model = "pooling")

non-spatial individual fixed effects (within) model using plm
FEwithin_nosp_crime <- plm(tot_crime, data = final_long_panel,
                  model = "within", effect = "individual")

non-spatial individual random effects model using plm
 REgls_nosp_crime <- plm(tot_crime, data = final_long_panel,
                model = "random", effect = "individual")

Hausman test using plm
hausman_nosp_crime <- phtest(tot_crime, data = final_long_panel)

Hausman tests robust to spatial autocorrelation using splm
hausman_sp_crime <- sphtest(tot_crime, data = final_long_panel,
     listw = queen, spatial.model = "error", method = "ML")

hausman_sp2_crime <- sphtest(tot_crime, data = final_long_panel,
        listw = queen, spatial.model = "lag", method = "ML")

So far so good. however, when i run the LM test for spatial lag dependence:
slmtest(tot_crime, data = final_long_panel, listw = queen,
        test = "lml", model = "within")

I get the following error:

Error in X %*% coef(mod) : non-conformable arguments.

Any idea what may be causing this error?
geoid           density     disad_inde  time    ndvi    rt_crime
170312909001    13035.08    2.33438992  1   0.31162748  97.70808203
170312909001    13035.08    2.33438992  2   0.49197907  106.1519903
170312909001    13035.08    2.33438992  3   0.4220205   100.1206273
170312909001    13035.08    2.33438992  4   0.13471775  95.29553679
170312909001    13035.08    2.33438992  5   0.3566384   95.29553679
170312909001    13035.08    2.33438992  6   0.4852591   97.70808203
170312909001    13035.08    2.33438992  7   0.48259588  90.47044632
170312909001    13035.08    2.33438992  8   0.17898722  92.88299156


texts: Is there a quicker way of checking if a feature exists with OGR

texts:
I'm writing lot's of features to a FileGDB using OGR.

The slow point in my script is checking that the feature isn't already in the Feature Class, based on the ID.  Currently I'm using a function and:

crfLayer.SetAttributeFilter("ID = " + str(crfID))
Count = crfLayer.GetFeatureCount()
if Count == 0:
    return 0
else:
    return 1

I've tried looping creating a list from the file and checking if the ID is in it, however this is slower.

Does anybody have any alternatives to the code above?


texts: await layer.LoadAsync(); invalid response | ArcGIS DevLab; Load an item from ArcGIS Online

texts:
So I am doing some excercises in the DevLab from the ArcGIS developer platform.
I did an excercise in which you load a layer from the ArcGIS Online portal. Their example code is the following (it works like a charm);

using System;
using System.ComponentModel;
using System.Runtime.CompilerServices;
using Esri.ArcGISRuntime.Mapping;
using Esri.ArcGISRuntime.Portal;

namespace LoadAnItemFromArcGISOnline
{
/// <summary>
/// Provides map data to an application
/// </summary>
public class MapViewModel : INotifyPropertyChanged
{
    public MapViewModel()
    {
        InitializeMap();
    }

    private Map _map = new Map(Basemap.CreateTopographic());

    /// <summary>
    /// Gets or sets the map
    /// </summary>
    public Map Map
    {
        get { return _map; }
        set { _map = value; OnPropertyChanged(); }
    }

    /// <summary>
    /// Raises the <see cref="MapViewModel.PropertyChanged" /> event
    /// </summary>
    /// <param name="propertyName">The name of the property that has changed</param>
    protected void OnPropertyChanged([CallerMemberName] string propertyName = null)
    {
        var propertyChangedHandler = PropertyChanged;
        if (propertyChangedHandler != null)
            propertyChangedHandler(this, new PropertyChangedEventArgs(propertyName));
    }

    public event PropertyChangedEventHandler PropertyChanged;

    private async void InitializeMap()
    {
        // Create a portal
        ArcGISPortal portal = await ArcGISPortal.CreateAsync();

        // Hold a portal item ID
        String portalItemId = "883cedb8c9fe4524b64d47666ed234a7";

        // Create the portal item
        PortalItem portalItem = await PortalItem.CreateAsync(portal, portalItemId);

        // Create the layer
        FeatureLayer layer = new FeatureLayer(portalItem, 0);

        // Ensure the layer is loaded
        await layer.LoadAsync();

        // Create the map with updated basemap and initial viewpoint
        this.Map = new Map(Basemap.CreateTopographic()) { InitialViewpoint = new Viewpoint(layer.FullExtent) };

        // Add the layer to the map
        this.Map.OperationalLayers.Add(layer);
    }

}
}

However, I wanted to try using a different layer. The layer I tried it with is a public layer created by the dutch government. It is a layer showing the heights of the netherlands. http://www.arcgis.com/home/webmap/viewer.html?webmap=4feca304f23847feae20eec8fec76fd8

So I simply replaced the ID of the tutorial with the ID from the Dutch layer and it gives me an error with the await layer.LoadAsync();

It sais the following:
Esri.ArcGISRuntime.ArcGISRuntimeException: 'Invalid response'

It does show a map with an esri basemap with a global extent. but the app is none responsive.

The entire excercise is done in Visual studio with ArcGIS runtime sdk for .net extension.

So what I'm thinking is that the dutch layer has a different coordinate system or maybe the layer has a maximum scale, or maybe the issue is that the data is raster. Does anyone has a tip/ know how to fix it, so that I can display this layer and maybe eventually more dutch layers?


texts: Using ArcGIS for clustering and correlating large datasets

texts:
I have a large amount of historical lightning strike data (with lat, lon and timestamp) and am trying to separate out isolated events, which aren't part of larger clusters found in thunderstorms (i.e. potentially anomalous events). I then want to see if the density/frequency of these isolated events is correlated to areas of earthquake activity (I also have a list of lat lon and timestamps for earthquake events.).

I have been trying to do this in python, but have been put on to the fact that GIS might be much better suited for the task. Could anyone tell me if this is so, and how I might go about it in a memory efficient manner?

I have access to ArcGIS 10.2.2 for Desktop.

An example of the sort of isolated event cluster I'm looking at can be seen in the bottom left of this image. As you can see it's in a smaller group, not attached to a larger storm. Sometimes there are only 1 or 2 strikes. my initial parameters for 'closeness' will probably be 200-300km, and 2-3 hours


texts: Converting .OBJ to .CTM file using openctm

texts:
I would like to convert a 3D .OBJ (Wavefront OBJ file) file to .CTM (OpenCTM 3D mesh file) format.
I know that openctm API (http://openctm.sourceforge.net/) have both Python and JavaScript compatibility. Here, I tried using openctm python API (https://github.com/lejafar/Python-OpenCTM/blob/master/README.md).
When I import .OBJ file, I am getting an error
OSError: Error loading file: b'CTM_FILE_ERROR'

When I try to import a .CTM file, it is working fine.
import openctm

# read
mesh = openctm.import_mesh('foo.ctm')

print(mesh.vertices.shape)
# (35947, 3)

# write
openctm.export_mesh(mesh, 'bar.obj')

So, I can import a .CTM file and export it as a .OBJ file. But, I am unable to import a .OBJ file.
What I want to do is to import a .OBJ file and export it as a .CTM file. How can I do this ?


texts:
I tried to compile postgis-2.1.1 in RedHat 6.3 with:

# ./configure --with-pgconfig=/usr/local/pgsql/bin/pg_config --with-gdalconfig=/usr/local/bin/gdal-config --with-geosconfig=/usr/local/bin/geos-config --with-xml2config=/usr/bin/xml2-config --with-gui  --with-raster --with-topology

It gives:

  PostGIS is now configured for x86_64-unknown-linux-gnu

 -------------- Compiler Info -------------
  C compiler:           gcc -g -O2
  C++ compiler:         g++ -g -O2
  SQL preprocessor:     /usr/bin/cpp -w -traditional-cpp -P

 -------------- Dependencies --------------
  GEOS config:          /usr/local/bin/geos-config
  GEOS version:         3.4.2
  GDAL config:          /usr/local/bin/gdal-config
  GDAL version:         1.10.1
  PostgreSQL config:    /usr/local/pgsql/bin/pg_config
  PostgreSQL version:   PostgreSQL 9.3.1
  PROJ4 version:        48
  Libxml2 config:       /usr/bin/xml2-config
  Libxml2 version:      2.6.26
  JSON-C support:       yes
  PostGIS debug level:  0
  Perl:                 /usr/bin/perl

 --------------- Extensions ---------------
  PostGIS Raster:       enabled
  PostGIS Topology:     enabled
  SFCGAL support:       disabled

 -------- Documentation Generation --------
  xsltproc:             /usr/bin/xsltproc
  xsl style sheets:
  dblatex:
  convert:              /usr/bin/convert
  mathml2.dtd:          http://www.w3.org/Math/DTD/mathml2/mathml2.dtd

But, when I compile I get:

# make
[...]
usr/local/lib/libgdal.so -L/usr/local/pgsql/lib -L/usr/kerberos/lib64 -ljpeg -lpng -lpq -lpthread -lrt -lcurl -ldl -lgssapi_krb5 -lkrb5 -lk5crypto -lcom_err -lidn -lssl -lcrypto -lxml2 -lz /usr/local/lib/libgeos_c.so /usr/local/lib/libgeos.so -lc -lm
../../liblwgeom/.libs/liblwgeom.so: undefined reference to `isfinite'
collect2: ld returned 1 exit status
make[2]: *** [raster2pgsql] Error 1
make[2]: se sale del directorio `/root/Descargas/postgis-2.1.1/raster/loader'
make[1]: *** [rtloader] Error 2
make[1]: se sale del directorio `/root/Descargas/postgis-2.1.1/raster'
make: *** [all] Error 1

Any idea?


texts: How to run Mapbox Vector Tile Mini Server?

texts:
I am trying to make an offline map for the vector tiles that I've exported from Mapbox Studio (mbtiles). I found this :

https://github.com/BergWerkGIS/mini-server

I have it downloaded on my desktop and it works online, but I can't get it to work offline. I just want to be able to have my vector tiles(mbtiles) to be displayed on the browser offline.


texts:
I'm new to GIS and I have 'part of' a problem to solve. Basically I have road details in geojson (WGS84) with LineString of points covering a 100KM square.
    I'm designing an architecture to snap a GPS co-ordinate to a grid cell (with different precision). So I have to construct a hash function which will generate a cellId (say 100m/1000m/10KM square) precision.

Say, if my X,Y is 58.345876, 10.765433 then applying this hash function should give me a cellId (4 bytes/8 bytes) with which I can look up. I tried different ways (likes of doing modulus etc.,) but it does not seem to be consistent.

Can someone guide me with code as to how to build this algorithm?


texts: Dynamically adding controls onto MapBasic dialog?

texts:
I have seen a lot examples with dialog creation, like:

    Dialog
    Title "Title"
    Width 198 Height 52
   Control Button
        Title "Обзор"
        Width 50 Height 15
        Position 133, 6
        ID 5001
        Calling SetTextEditFolder
    Control Button
        Title "Создать"
        Position 15, 28
        Width 50 Height 15
        Calling CallforYes

But is just a static dialog.  It's possible to dinnamically add anf remove controls from dialog?  For example i want press button or select value in list edit and change set of controls in dialog.
Any advises?


texts:
I understand that the result of InSAR interferometry yields a LoS deformation and if deformation in the vertical or horizontal is to be computed, it should require additional gps observations at the site.

My first question is, could I simply solve for the components (vertical and horizontal) using geometry from the incidence angle?

Secondly, for DEM generated by repeat pass SAR images, i understand that the height obtain is relative. Can someone explain to me what does relative mean? Does it refer to the change in height from the first image to that of the second image? How then can I get an absolute height?


texts: Problem starting QGIS 2.14, Prozedureinsprungpunkt is not found

texts: I have a servere problem concerning starting QGIS 2.14.3. I haved installed that version using OSGEO4W installer. Besides, I have QGIS 2.8.6 on my machine to keep the old LTR version. As it is mentioned in the error box the DLL could not find the right version, isn´t it? I checked for duplicates on dll´s in the windows system folders (even the hidden folders) - no success! Does anybody has an idea how I could solve that problem?


texts: Progress dialog box of ArcGIS Python script tool stealing focus from TKinter GUI?

texts:
I created a script tool in the custom toolbox, it runs like this:

But the default progress dialog keep stealing focus from the Tkinter GUI, that makes me so hard to move/resize the OpenFileDialog, expand the dropdown menu, double-click to choose item... ...

I have to right-click on the file to select it, paste the path into the inputbox to change the location.

Is there any way to avoid this or disable the default progress dialog?


texts: QGIS 1.8 startup crash - looking for a log file

texts:
I have a problem with QGIS 1.8 with OSGEO4W installer (windows 7 - 64 bit).
QGIS crash on startup, it shows the splash screen and then nothing.
I have updated from an old 1.7.4 QGIS version.
I have done the same procedure in a portable computer and it works fine.
The problem is only in the workstation.
Is there a log file to view what is the problem?
Thank you


texts: I want to arrange 96 points as a circle/ring. Is there a way to make this automatically, for example "Buffer", "vertices"?


texts:
I have geojson data that draws shapes similar to the following:

I wanted to know if there are any operations that can give me the geojson which represents the outline of these shapes represented by the red lines in the following image:

The coordinates do not have to match exactly what is represented by the red lines I only need the coordinates to wrap the outside boundaries of the shape.


texts: How to remove or reset all metadata from FileGeoDatabase (*.gdb)?

texts:
Removing metadata from shapefiles is fairly easy, since they are stored in a seperate xml file.

How can I remove or reset all metadata in a FileGeoDatabase (*.gdb)?


texts: ST_Overlaps in return TRUE for adjacent polygons - PostGIS


texts: PostGIS raster geoTiff simple Lat/Lon query example

texts:
I've successfully loaded a geoTif into PostGIS (Postgresql 9.3.9) using the following command on my Ubuntu 14.04 server:-

sudo -u postgres raster2pgsql -s 4326 -d -I -C -M -R -l 4 /path/to/my_raster.tif -F -t 10000x10000 my_raster | sudo -u postgres psql -d my_database

It's a 1 banded grayscale raster (lat/lon) with cell values representing some environmental factor on the ground. Hence, I simply want to get the cell value at a given lat/lon.

Following some examples I found I came up with the following SQL:-

SELECT
  ST_Value( rast, ST_SetSRID(ST_MakePoint(0,51.5),4326) ) AS rast_val
FROM
  my_raster
WHERE
  ST_Value( rast, ST_SetSRID(ST_MakePoint(0,51.5),4326) ) IS NOT NULL;

This seems to get the right value for me when I cross reference the same location loading the raster in QGIS.

My question is basically a sanity check if the above SQL is correct as I am not 100% sure of the WHERE clause.


texts: Measure Line while editing in Edit Mode (Workaround?)

texts: It's unfortunate that QGIS does not allow the use of the Measure Line tool while activiely editing something with Add feature or Vertex tool. As soon as one of the latter tools are switched on the measured line disappears along with the measure window. In some cases however it is necessary to place a point or a line in a certain distance off another object. Therefore it would be handy to see the measured line while editing. Is there a workaround maybe or an alternative to the standard Measure Line tool or maybe any configuration which allows the Measure Line tool to stay open during editing?


texts: How can I include metadata for my data in PostGIS? I want to have the associated metadata for my shapefiles saved alongside the data in my database so that I can connect the database to QGIS and display the data there.


texts: How to read shapefile from .zip with OGR command line tools?

texts:
I have thousands of shapefiles in .zip archives to work with. How can I use ogr2ogr and other GDAL/OGR tools to work with these data without having to extract the archives first?

From gdal cheatsheet and Reading a GDAL dataset in a .gz file or a .zip archive I know that I can use the visizip driver, but for formats that consist of multiple files per data layer like shapefile you need to explicitly name the file inside the archive:

ogrinfo /vsizip/archive.zip/105A_wetland_a.shp

How might I efficiently work with archives, that could contain multiple shapefiles per zip, without having to know the and pass the name of each layer ahead of time? e.g.

ogrinfo /vsizip/archive.zip

Example archive file: bndt_105a_shp_en.zip from Geogratis.


texts: Seeking GIS that can calculate fractal dimensions?

texts:
As ‘Fractals are Everywhere’ in geo-ref data, which GIS other than ArcGIS
can calculate fractal dimensions?


texts: What is the best way to create a line of sight simulation raster?

texts:
I have been tasked to create some sight-line simulations for a powerline project. I have several (85) observation points and need to create a simulation of what the view would look like, looking from the observation points back to a proposed powerline.

I have a shapefile x,y,z of observation points, photos from observation points to proposed powerline, and a 3d drawing of powerline in autocad.

Any suggestions on the best work flow?

I have a ArcGIS licence with 3d extention, and Autocad Civil 3D.


texts:
How to combine below expressions

Function FindLabel ( [Medium], [LocationStatus], [Status] )
  if  isnull (  [Medium] ) then
    FindLabel= ""
  else
    FindLabel = [Medium] & " Flowline"
  end if
End Function

Function FindLabel ( [Medium], [LocationStatus], [Status] )
  if  isnull  (  [LocationStatus] ) then
    FindLabel=""
  else
    FindLabel = "    Location-" & [LocationStatus]
  end if
End Function

Function FindLabel ( [Medium], [LocationStatus], [Status] )
  if  isnull (  [Status] ) then
    FindLabel=""
  else
    FindLabel = "    Status-" & [Status]
  end if
End Function


texts: Error while comparing two raster layer in R

texts:
I have got a geographic point distribution with N prehistoric settlements with a size of x and want to compare the distribution with an economic simulation, which produces M settlements (Usally N > M). I've done about 110 runs with different parameters and now I'm looking for a best-fit model. Therefore I calculate the Kernel Density Estimation for the real settlemts (x_dens), the simulation results (y_dens) and for the M-biggest real settlements (xM_dens) and want to compare the cell values of two kernel density estimations in R.

As a first step, i am using the function density.ppp from the Package "spatstat".

x_dens <- density(x_ppp, sd, eps=rw, edge=TRUE, at="pixels")
x_ras <- raster(x_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

y_dens <- density(y_ppp, sd, eps=rw, edge=TRUE, at="pixels")
y_ras <- raster(y_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

xM_dens <- density(xM_ppp, sd, eps=rw, edge=TRUE, at="pixels")
xM_ras <- raster(xM_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

The next step includes the calculation of the cell differences for x with xM, to see which would be minimum error for my M simulated settlements:

xxM_comp <- sqrt(((x_ras - xM_ras)^2 / x_ras))
xxM_comp [!is.finite(xxM_comp )] <- 0
xxM_comp_result <- sum(xxM_comp @data@values)

As a next step, i'm doing the same for my Y-simulated settlements:

xy_comp <- sqrt(((x_ras - y_ras)^2 / x_ras))
xy_comp [!is.finite(xy_comp )] <- 0
xy_comp_result <- sum(xy_comp @data@values)

To evaluate the best-fit model, i caluclate:

results <- xy_comp_results/xxM_comp_result

xy_comp shall result in a new raster layer, which shows the differences in both point pattern, but the resulting cell values are enormously high at some points and as a consequence xy_comp_result has a illogically high value. It is due to the fact, that some pixel values of x_ras much to small compared with y_ras. Do you have got any suggestions, how to solve this problem. I've tried to remove outliers, but I did not get proper results.

The methodology is described in D. Stelder, Where do Cities Form? A Geographical Agglomeration Model for Europa, Journal of Regional Science 45 (4), 2005, 657-679 (http://www.regroningen.nl/stelder/doc/JRS_nov2005_c.pdf) Starting with page 669.


texts: Making buffer from line using GDAL and Python

texts:
I am trying to make a buffer from a line using gdal.

I use the example in the gdal Cookbook

import ogr, os

def createBuffer(inputfn, outputBufferfn, bufferDist):
    inputds = ogr.Open(inputfn)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputBufferfn):
        shpdriver.DeleteDataSource(outputBufferfn)
    outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
    bufferlyr = outputBufferds.CreateLayer(outputBufferfn, geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()

    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)

        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)

    def main(inputfn, outputBufferfn, bufferDist):
        createBuffer(inputfn, outputBufferfn, bufferDist)

    if __name__ == "__main__":
        inputfn = 'obli.shp'
        outputBufferfn = 'oblibuffer001.shp'
        bufferDist = 0.01

        main(inputfn, outputBufferfn, bufferDist)

 # creating the projection file. Not included in the cookbook
    prj = open('oblibuffer001.prj', 'w')
    proyeccion = 'GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],' \
             'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
    prj.write(proyeccion)
    prj.close()

It works fine

http://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html?highlight=buffer

But I can not understand exactly what kind of units I must use in the bufferDist

Ussing 0.01 I think it was a kilometre but the result is different depending the direction of the line. For example
Using a parallel_direction_line using  bufferDist = 0.01 meassured in the map the buffer is 1,1 km
       a meridiam_direction_line  using  bufferDist = 0.01 meassured in the map the buffer is 0.456 km

but  a oblicual_direction_line  using  bufferDist = 0.01 meassured in the map the buffer is 0,57 km

How can I find the value for the bufferDist variable I think depends on the latitude and the direction of the line?


texts: How to enable Python scripts in processing toolbox in QGIS 3.16

texts:
I'm not able to enable Python scripts in the processing toolbox. It should turn up in the toolbox (see the third screenshot).
The Python logo is also missing from the toolbox.
Is there something wrong with the settings in the second screenshot?
I tried to write my own custom script. When that failed I exported a model from the models as 'Python Algorithm' but also no success.


texts:
I am tring to use OSM2PO on Windows 7 but nothing works.
I have JAVA Already installed. When I launch the command
java -jar osm2po-core-4.7.7-signed.jar prefix=at "USTHB.osm"

I get this :

AnySolutions ?


texts: QGIS Advanced digitizing switch from Cartesian to Ellipse

texts:
I´m trying to draw a line per length, but I need to enter the length in the Ellipse value.

Is there a way to switch the values in the length window in the advanced digitizing panel from Catesian to Ellipsoid?


texts: How do I combine multiple raster layers to make a landcover map for Borneo?

texts:
I am currently making a landcover change map for a Kinabatangan, Sabah, Borneo for my undergrad dissertation. I downloaded landsat MSS and 8 images. These were pre-processed and both classified using unsupervised isodata classifications. Due to the restricted resolution I had to then use Google Earth Pro to delineate and map palm oil concessions (large palm oil sites) and smallholders (small farmers). I made polygons for each type and converted these to KML for GIS. These polygons were then changed to raster format so everything was the same. Both rasters were combined to make one palm oil layer. I standardised all the data e.g. all had the same projected coordinate system and cell size.

What I can't figure out though is how to combine my palm oil raster data with my 2014 landcover map (this was the only year the palm oil was seen). I essentially want three maps

1) 1973 Base landcover map, with these classes: Water; Forest and Non-forest

2) 2014 landcover map with combined palm oil data, with these classes: Water; Forest; Non-forest; Smallholders; Large Concessions

3) 1973 - 2014 Landcover change map, with these classes: Water; Forest; Non-forest; Smallholders; Large Concessions.

I have tried to use the raster calculator and also merge to new raster tools. Any comments of what I can do / am already doing wrong would be appreciated! P.S. I have attached some pictures of what I am working on [1: Overview of map; 2: Overview of classes].


texts: OpenLayers 4 - Preserve geometry shape when modifying

texts:
Using the ol.Draw interaction and ol.interaction.Draw.createBox() geometry function we can draw a regular rectangle:

let geomFunction = ol.interaction.Draw.createBox();
let draw = new ol.interaction.Draw({
  source: layer.getSource(),
  type: 'Circle',
  geometryFunction: geomFunction
});

but when using ol.Modify interaction to modify the rectangle I can't ensure the geometry will still a regular rectangle because I can't set a geometry function to the modifying interaction!
How can I ensure or preserve the geometry shape when modifying?


texts: Naming raster layers in a rasterstack with the same name

texts:
I have a for loop where I loop layers from a rasterstack into a predict function, which requires the output of a model and two raster files, including one layer for the present ("rastFile") and one stack or brick for the future ("stackedRast"), as follows:
  for (i in 1:94){

  timePred <- predict(G15.gdm.1, rastFile, time=TRUE, predRasts=stackedRast[[i]])

  stack(timePred,stackedRast[[i]])

  writeRaster(timePred, "xxxxx.tif")

  rm(timePred)

}

First of all, for the predict function to work, I need all the layers in the rasterstack to have the same name, for instance "bio1"
for (i in 1:nlayers(stacked)) {names(stacked[[i]]) <- "bio1"}

However, R adds a number to the end of the layer names, like this:
> stackedRast
class      : RasterStack
dimensions : 1200, 2760, 3312000, 4  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -180, -65, 35, 85  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names      :    bio1.1,    bio1.2,    bio1.3,    bio1.4
min values : -24.70417, -22.44583, -19.79583, -16.23750
max values :  25.96667,  27.47917,  28.76250,  30.28750

Is there a way to get around this? (P.S. I've just included a smaller rasterstack of 4 layers for the example, but my final rasterstack will have 94 layers).
Does anyone please have a suggestion for how I can elegantly create one process (loop or otherwise) for  changing the layer names so that they all match?


texts:
I'm trying to run GeoPackage on QGIS Server and query the information. I created a file through GeoPackage's QGIS Desktop.
To do the queries I used the following tutorial:
https://docs.qgis.org/3.4/en/docs/user_manual/working_with_ogc/server/getting_started.html#installation-on-debian-based-systems
I am trying to use Nginx and did not understand exactly where in the configuration file I am adding the next part?

Also, how do I properly link GeoPackage to a Server?


texts:
I need to combine 2 large polygon datasets. One of these datasets contains more recent data, so the polygons from this dataset need to be copied into the existing dataset.

There are three different instances where the existing dataset would need to be modified. See below for screenshots and explanations.

(In all screenshots the dashed red outline represents the existing data. The green polygons are new data.)

Situation 1 - Entirely new polygon:

Action - New polygon needs to be added to existing dataset. Polygons that remain the same can be ignored.

Situation 2 - An existing polygon has changed shape:

Action - Existing polygon(s) needs to be deleted. New polygon then needs to be added to existing dataset.

Situation 3 - Multiple new polygons within an existing polygon

Action - Existing containing/large polygon needs to be deleted. New polygon then needs to be added to existing dataset.

So, my question is; would it be possible to combine these datasets according to the rules above, either in QGIS or PostGIS?


texts: Folium, make legend disappear when switching layer (LayerControl)

texts:
I created some legends in a map using Folium and Python using the html example by InLaw(1) in this answer. From that answer, I also manage to get this legend using a MacroElement from Branca (2), where template is an html macro template. What I want is that any of those legends appear and disappear from the map in accordance to a LayerControl.

(1) map.get_root().html.add_child(folium.Element( legend_html ))

(2) macro = MacroElement()
macro._template = Template(template)
m.get_root().add_child(macro)

What have I tried? I tried adding the html and MacroElement to the FeatureLayer or a Choropleth map with no avail. After lurking for a while in the web I found this example that uses java script (I guess, I have -0 experience with JS) to make the bar invisible. But, when I substituted the LinearColormap for any of the above, the legend did not disappear. For example, using the variable macro from (2) I tried this:

m.add_child(BindColormap(FeatureGroup, macro))
m.add_child(BindColormap(Choropleth, macro))

Where FeatureGroup and Choropleth are the names of the variables named like that for simplicity for this example.

I don't know what else I can do.


texts: Converting symbol style set for ArcGIS Desktop for use in QGIS?

texts:
Is there a way to convert a symbol style set for ArcGIS for use in QGIS?

There are symbols used in Incident Command I would like to use in QGIS, but can only find a style set for use in ArcGIS.  The style set I have is saved as a .style file.


texts: Why can't I "Add label columns to attribute table" etc. in QGIS 1.9?

texts: Why can't I "Add label columns to attribute table" and why can't I use "Letter spacing", "Word spacing" and so on? I have a shp-file. QGIS version 1.9


texts: Field calculator for output line from ModelBuilder?

texts:
Calculate values for only a single feature in Model Builder output.

How do I calculate a field for only the newly generated feature? In my  model everything goes fine but the calculate field operates on all the features in output layer not just the newly created feature.

I am attempting to achieve this using a selection in my model.


texts: Faster way than Spatial Join to join field from attribute table to address table?

texts: I have a table of addresses for every state in the nation. I just recently finished doing spatial joins and adding in census block/tract data per state. I have a national shapefile of CoreBasedStatisticalAreas, is there a faster way to add one field from this attribute table to my completed address tables (millions of records) than doing Spatial Joins?


texts: I want to understand the concept of a "value raster" and if I can use it in neighbourhood statistics.


texts: I need to implement Access Control List (ACL) for my web application and I need to know common practice of using permissions in GIS systems. What levels of acces exists: layer, feature, attribute, extent, maybe anything else? What types of activities can be used by users: editing, browsing? Any response would be appreciated.


texts: Problem with dll while installing QGIS in windows

texts:
I am having an issue installing QGIS on my windows 10 machine. It used to be finely installed, but due to some reasons I had to uninstall.
It has to do with older (2.16) and newer versions (3.0,3.10..).
I am getting a dll related error as you can see in the image attached.


texts: How save an layer styles as Default in ArcGIS for Desktop?

texts: is there any way to save the layer styles as default in ArcMap?


texts: Distance with Bearing Calculation not working

texts:
I have this function in Python that calculate distance from a point with bearing. Problem is the latitude part is generating the wrong answer and i cannot figure out where i went wrong:

def add_distance(self, lat, lon, bearing, distance):

    d = distance/6371
    x = sin(lat) * cos(d) + cos(lat) * sin(d) * cos(bearing)

    new_lat = asin( float(x) )
    new_lon = lon + atan2(sin(bearing)*sin(distance/6371)*cos(lat),
                     cos(distance/6371)-sin(lat)*sin(new_lat))
    print '\tNew Lat:%s' % (new_lat )
    print '\tNew Lon:%s' % new_lon
    return ''


texts: What are convenient uml tools to create geodatabase models?

texts:
I used to model databases with Microsoft Office Visio, then export the model to xml and import geodatabase schema to GDB with ArcGIS tools.
Since Microsoft no longer support Visio export to xml, what are the alternatives?

Has anyone used Enterprise Architect? Is it convenient for designing Geodatabases for ArcGIS?

I have tried ArcGIS Diagrammer (for ArcGIS 9.2/9.3 or for ArcGIS 10) but it is not so convenient..


texts: Extracting values from Corine Land Cover (CLC) 2018 raster using R

texts:
I would like to extract values from Corine Land Cover (CLC) 2018 raster.
library(sp)
library(rgdal)
library(raster)

r <- raster('U2018_CLC2018_V2020_20u1.tif')
#Warning message:
#In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#  Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition

I got a Warning message, probably due to a lack of a +datum parameter in the proj4 string.
The Coordinate Reference System (CRS) for the CLC 2018 is EPSG:3035 (ETRS89, LAEA).
The crs from the raster 'r' is :
crs(r)
#CRS arguments:
# +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

Then I load my 'data' file that is a 'data.frame' with 419389 rows and 15 columns including the 'Latitude' and 'Longitude'.
I then run these few lines
lats <- c(data$Latitude)
lons <- c(data$Longitude)
coords <- data.frame(x=lons,y=lats)
points <- SpatialPoints(coords, proj4string = r@crs)
points <- na.omit(object = points)
values <- extract(r, points)

points
#class       : SpatialPoints
#features    : 419389
#extent      : -8.06205, 2.39726, 49.18919, 60.59815  (xmin, xmax, ymin, ymax)
#crs         : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

However, the 'values' return NAs.
I guess the problem is the projection, so I tried to reproject the raster but it doesn't work.
projection(r) <- crs("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +datum=WGS84 +units=m +no_defs")


texts: Layer use too big for the coordinate in GE

texts:
I've got a problem here. I already changed the format from Tab file into Shp file using MapInfo suggested by friends here.

My question is, after we set all the layer into the right coordinates, even using Google Earth to set the coordinates, but when added as a new layer, the size of the layer added seems bigger than the one in the google earth layer (meaning that the Google Earth layer is so small when it extent)

How can I change the size of the layer so that the point in Google Earth matches the one with the layer added afterward?  Do I need to re-size the layer or another way?  Someone please help!
Thank you...

1

My 1st layer capture from GE using mapwindow GIS

2

I didn't click the full extend button yet. This is my 2nd layer and it shows the 2nd layer was purple in color were actually a mega size layer. i think around more than 50 times bigger than the first layer.

3

This layer after full extend. and i cannot manage to find the 1st layer at all because it is too small.

Q: How can I change the size of the 2nd layer so that it will be the same size as the first layer with the same coordinate?


texts: Extract Values to Point from Raster Issue

texts:
I am trying to use the tool Extract Values to Points on a File Geodatabase raster, but I am getting some weird errors.

I've used this tool multiple times before with no issue, but not since I upgraded to 10.1, so I'm wondering if there is an issue with it.

The errors I got back so far have been:

This error happened when I specified a shapefile as the output: The specified geodataset is invalid for the conversion. The syntax help for the tool you are using should indicate the supported types of input data. For example, for the particular conversion, a rasterlayer may be an invalid type of geodataset.
Next I have an error about it needed at least one band in the raster, the raster does has one band. This occurred when I specified a File Geodatabase as the output.
Finally I had the good ole Error: 999999 for no apparent reason.

The raster is a DEM derived from LiDAR that was not created by a different company. It's too large to be made available for testing, but here are its specs:

Raster_Information

ColumnsandRows  5119,8241
NumberofBands   1
CellSizeXY      1, 1
UncompressedSize    160.93 MB
Format      FGDBR
Source_Type     Generic
Pixel_Type  floating point
Pixel_Depth     32 Bit
Colormap    absent


texts: ArcObjects- IBorder.GetGeometry() in Python

texts:
I am trying to get geometry of frame border trough IFrameElement.Border. I do receive IBorder object, which has .GetGeometry() method.

How do I use this method to get geometry of the border?

I can not find code snippet in any language, so any advice where to look for is appreciated.


texts: Counting number of unmasked pixels per image in collection using Google Earth Engine

texts:
I'm trying to count the number of pixels that are masked by the cloud mask of the image collection given by the Earth Engine dataset. I have clipped my collection with a ROI and now I want to keep the images with 10% or less of clouds  (masked pixels). By now, I am only able to count them for a single image but not for the entire collection.

I am now using this code for one image but I am not able to reproduce this for the entire collection.

var Peyto_pixelscount = Peytoclip.first().select('ALBEDO').reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: Peyto_glacier.geometry(),
  scale: 30,
  maxPixels: 1e9
});

print(Peyto_pixelscount)


texts:
I'm using owslib to download an aerial photograph. In the WMS GetMap request I specify a srs, so I assume (maybe wrongly?) that the response I get contains a crs.

When writing the result to a GeoTiff the file doesn't contain spatial information. Is there a way to do this or do I have to use GDAL to project the result, and if so, how?

This is the code I'm using so far:

from owslib.wms import WebMapService
import gdal
from osgeo import osr

wms_url = 'https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wms?'
wms = WebMapService(wms_url,version = '1.3.0')

img = wms.getmap(layers=['Actueel_ortho25'],
                 srs='EPSG:28992',
                 bbox=(97712,440638,98220,441091),
                 size=(1000,1000),
                 format='image/GeoTIFF',
                 transparent=False)

out = open('path\\to\\test.tiff','wb')
out.write(img.read())
out.close()

print('done')


texts: Pasting coordinates into QGIS 3.x from clipboard

texts:
I was looking for a case to easily add points to a layer from coordinates. I my case I often want to quickly add one or a few points manually to an existing layer. I have seen an occasional reference that it would be possible to paste it directly into a layer. It would fit perfectly in my workflow to copy some coordinates in an external application e.g. google maps/osm to paste (ctrl-v) into an existing layer.
So I tried pasting the syntax in the two CRS used in my current file:

"point(571287,6799222)"  (CRS: New Amersfoort)
"point(48.0578195,7.2497044)" (CRS: WGS84)

Than QGIS prompts that features are pasted succesfully:

The problem is that indeed a extra row is added to the table, but no coordinates are actually added. Also when I 'extract' the coordinates in the table field calculator with the functions $x adn $y it appears as null instead of points added in a different way and which are showing or giving an output :

I guess there is something wrong in the syntax but I coul not find any references or guides on what syntax QGIS is expecting. Is there somebody who can provide more background on how this features is supposed to work?

P.s. I know there are different ways of getting points of exact coordinates into Qgis. My question is that apparantly this low-key workflow is supposed to work and I cannot figure out how or found any more info on this.


texts: Add a canvas on the map that can pan and zoom along with it

texts:
I want to know if it's possible to add a canvas on a ArcGIS map that can pan and zoom along with it.

The reason I ask that it's because we have different workflows where we have a lot of entities (+1000) that need to be refreshed at a regular interval (3-5 sec). These entities are SVG DOM element and the rendering performance is poor.

One solution we want to try to improve rendering performance is to add canvas on map and add these entities on it.

Is it possible? Do you have any idea/sample how to do that?


texts: i.landsat.toar doesn't Work With Band 6 of Landsat 7 in Grass 7

texts:
I'm facing an issue when running i.landsat.toar with a Landsat 7 image, here's the list of files:

LE70010752014040CUB00_B1
LE70010752014040CUB00_B2
LE70010752014040CUB00_B3
LE70010752014040CUB00_B4
LE70010752014040CUB00_B5
LE70010752014040CUB00_B6_VCID_1
LE70010752014040CUB00_B6_VCID_2
LE70010752014040CUB00_B7
LE70010752014040CUB00_B8

When the command is processing the bands, it gets stuck in band 6 and gives this error message:

ERROR: Unable to open header file for raster map <LE70010752014040CUB00_B61@>.

The rest of the band (7 and 8) aren't processed, neither with reflectance or radiance selected.

I configured the input_prefix=LE70010752014040CUB00_B and output_prefix=reflectance_


texts: How to solve error: 'user-defined methods must return a value' in GEE?

texts:
I would like to iterate a map function over several list-of-lists.
Here is the code I run:
var feature1 = function(h) {
    var FCol= ee.FeatureCollection(ee.List(losslist.get(h)).map(function(element){
        element = ee.Dictionary(element); // cast
         var list = ee.List(element.get('groups'));
         var feats = ee.FeatureCollection(list.map(function(listEl){
         var dict = ee.Dictionary(listEl);
     return ee.Feature(null, dict);
  }));
   return feats
}))};

var myList = ee.List([lossperelevation,lossperslope,lossperyear]);
var listNumbers = ee.List.sequence(0,2)
var result = listNumbers.map(feature1);

As an output, I get an error: User-defined methods must return a value. What I understood is that I need to add another return line. I do not understand what it should return, as the output I need is feats.


texts: Is there a way to paste geometries in a workspace at the location i want? e.g. i want to paste a circle from another workspace in the workspace i'm working with. I need only the geometry, not the coordinates.


texts:
In this question I was asking for a aggregate solution to update a attribute (ID) if I change or add an new feature and got a great answer from OGmaps:
Add and change attribute / feature cascaded (QGIS)
This is simple and perfectly working! But I found out now, that this is consuming too much resources if I have more then 300 polygons. So I need a other way and thought about SQL - I'm sure this can do the same job maybe with more performance. So how do I translate this aggregate function:
aggregate(
  layer:= 'Polygons',
  aggregate:= 'concatenate',
  expression:= Id,
  concatenator:= '-',
  filter:= contains($geometry, geometry(@parent)),
  order_by:= -area($geometry)
)

Into a SQL-Expression?

Here I uploaded a example data set: https://www.map-site.de/tmp/sql-cascade.zip I need the green dynamic ID's (with place holder).

The green ID's are generated by connected aggregate functions in a virtual field. If any object name is changed, the ID will update - very nice but too slow for a good workflow (I expect up to 1500 polygons)


texts: Why can't I find the TIN Management tool in 3D Analyst toolbox?

texts: I am looking for the CREATE TIN tool in the 3D Analyst toolbox, but I can´t find it. Any ideas why I can't find it?


texts: Regarding Atmospheric correction and Elevation DEM Projection system

texts:
I am trying to conduct atmospheric adjustment of Landsat 8 image bands using 6s algorithm (i. atcorr grass gis tool). It seems, in the input section, an elevation raster map should be provided. In my literature review, I could not find any reference discussing the DEM's relevance in the process. It is not a required option either. Any explanation how and where it fits in the overall process? In my control file, I am already providing average height above msl (of surface I want to correct). My earlier attempts at correction led to max values above 1 (which is unusal, I was told). Want to explore if the DEM makes any difference or not. Anyways, I may need the DEM for topographic corrections later.

Correction: Earlier, While generating Radiance values from DN number using i.landsat.toar; Atmospheric correction method at Metadata section was set to one of the options. It should be left blank while generating Radiance. This lapse may have contribute to later 6s Algorithm adjusted bands to have max value higher than 1.

I am trying to be consistent with the data. The Elevation source in line 12 of metadata references GLS2000. Online query led me to GLSDEM. I am assuming this is the DEM used in the L1T correction.

The original DEM is in WGS1984 geographic coordinate system. However, the image bands are in WGS1984 UTM projected coordinate system. Do I need to convert the DEM from geographic to projected coordinate system at all? I am asking because conversion will require bilinear interpolation and this will introduce distortion in the elevation values.

Update: Original Geometric projection GLSDEM do not even coincide with TOAR processed Radiance Bands once imported in Grass GUI. Most likely, have to switch to projected coordinate system.


texts: Shifting in pixel after subsetting with gdalwarp

texts:
I am trying to subset multiple image to smaller Area of interest
I have tried with gdal_translate & gdalwarp in commandline but it causes slight shift in pixel in the output raster after cropping/sub-setting image into smaller place by using 1st command.

my 1st command is

gdal_translate -projwin 5.83 47.86 10.98 45.65 -of GTiff inputraster.tif outputraster.tif

I have also tried in qgis (which is basically same thing) by clip raster by extent which also causes shift in pixel from original.

My 2nd command

gdalwarp -tr 0.01 0.01 -te 5.83 47.86 10.98 45.65 -s_srs "EPSG:4326" -t_srs "EPSG:4326" -r near -tap -of GTiff input.tif output.tif outputraster.tif

which causes changes in pixel value the highest raster value decreases in 210 from 215 in original raster


texts: Send Feature Attributes Along with Geometries to Difference Operation of Geometry Service

texts:
I've large number of data(geometries and attributes). I've to find their difference with other geometry, and after difference is calculated I've to assign attributes to those geometries. Example

Attribute-A  Geometry-A
Attribute-B  Geometry-B
Attribute-C  Geometry-C

and

SubtractorGeometry-B

[Note: Geometry-B will get subtracted by SubtractorGeometry-B]

Expected Result:

Attribute-A  Geometry-A
Attribute-C  Geometry-C

The problem I'm facing is that difference operation only accepts geometries and not attributes. When difference returns result, I don't know which attribute do I've to assign to which geometry.

One workaround is to send each geometry separately and assign attributes when it returns result, but for large data that's impractical  as for 1000 records their will be 1000 service calls.

Can anyone please suggest any solution?


texts:
Where can I find a comprehensive crosswalk for FIPS codes at the Place level?

I need to match Census Places from different data sets spanning the period 1980-2010, where each data set is using FIPS code as a Place identifier.


texts: problem in spline with barriers in ArcGIS 9.3?

texts: I am trying to use the script tool "Spline with Barriers" from arcgis 9.3 3D Analyst toolboxes. I am not getting expected output. The output is just only one value over the surface. We are using ArcGIS 9.3 in XP Remote Desktop and it is a German Version of Windows, so the system represents the decimal place as a comma. Is it possible to edit this tool to store the field values as 3.40 instead of 3,30?


texts: Export multiple Layer Definition Files (qlr)

texts:
I have a postgresql database with approximately 400 tables. I need to generate Layer definition files (qlr) in QGIS for each individual table. Making this manually would take huge amount of time, so I would like to make a batch process or a python script, that does it for me.

Is there anyone here that can help me in the right direction? :)


texts: Only using some data from shapefile in QGIS

texts:
How can I use only a few of the areas outlined in a shapefile? In this case, my shapefile has all of the ZIP codes for a U.S. state, and I would only like to display the ZIP codes for a particular city. Is there any way to do this?

Also, some of the ZIP codes are partially in the city, and partially out of it. How would I go about editing the boundaries of the outlined areas?


texts: Creation of tiled grid in JAVA with geotols

texts:
I am trying to create a grid using GeoTools in JAVA, covering specyfic area with resolution 50 m (distance between two points). To do so I have two points: the most nortwestern (let's call it NW) and the most southeastern one.

I created a table and put the NW point in tab[0][0]. Starting from tab[0][0] using the Geodetic Calculator (link to the doc) I am calculating the next points and putting them in right place into my table.

Here is the code I developed:

public void createGrid(double granularity) {
    logger.info(GridManager.class + " creation of grid called");

    int defaultCellID = 0;
    float defaultSigLevel = -1000;

    /*  __ __ __ __ __
     * |NW|__|__|__|NE|
     * |__|__|__|__|__|
     * |__|__|__|__|__|
     * |SW|__|__|__|SE|
     */

    // mind, this is lon/lat
    Point2D.Double nePoint = new Point2D.Double();
    nePoint.setLocation(sePoint.getX(), nwPoint.getY());

    Point2D.Double swPoint = new Point2D.Double();
    swPoint.setLocation(nwPoint.getX(), sePoint.getY());

    // calculations of a width of the grid
    this.gc.setStartingGeographicPoint(nwPoint);
    this.gc.setDestinationGeographicPoint(nePoint);
    double d = this.gc.getOrthodromicDistance();
    int width = (int) Math.ceil(d / granularity);

    // calculations of a height of the grid
    this.gc.setDestinationGeographicPoint(swPoint);
    d = this.gc.getOrthodromicDistance();
    int height = (int) Math.ceil(d / granularity);

    // grid creation
    grid = new TileCartesian[height][width];
    grid[0][0] = new TileCartesian((float) nwPoint.getY(), (float) nwPoint.getX(), defaultSigLevel, defaultCellID, 0, 0);
    double dist;

    for (int i = 0; i < height; i++) {
        this.gc.setStartingGeographicPoint(grid[i][0].getLon(),grid[i][0].getLat());
        dist = granularity;

        for (int j = 1; j < width; j++) {

            this.gc.setDirection(90 /* azimuth */, dist);

            Point2D.Double p = (Point2D.Double) this.gc
                    .getDestinationGeographicPoint();

            grid[i][j] = new TileCartesian((float) p.getY(),
                    (float) p.getX(), defaultSigLevel, defaultCellID, i, j);
            dist += granularity;
        }

        if (i != height-1) {
            this.gc.setDirection(180 /* azimuth */, granularity);
            Point2D.Double p = (Point2D.Double) this.gc
                    .getDestinationGeographicPoint();
            grid[i+1][0] = new TileCartesian((float) p.getY(), (float) p.getX(), defaultSigLevel, defaultCellID, i+1, 0);
        }
    }

    logger.info(GridManager.class + " Grid created");
}

My problem is that the most SE point of the calculated grid is far from the assumed one. And when I am exporting the result to GeoTIFF this cause deformation.

I really do not know how to solve this problem. Maybe there exists any other way to create such a grid, but I couldn't find any solution which would fit my needs.


texts: While doing the parallel line method by using the cadtools its asking the distance in meters but i need the parallel in feet values(I need to enter the distance in feet values). currently we are using the project file with the coordinate system as WGS84 EPSG:4326. I have changed the layer projection for feet values but its not working properly. can you please suggest on this that how can i change the units from meters to feet without changing the projection of project file.


texts: Why the shapefile of worldwide country boundaries is so large (>500 MB)?

texts: I downloaded GADM administrative level 0 worldwide and it seems very simple file: 255 rows of polygons representing countries. But the size of the file is very large: >500 MB. Why it's happening? Usually table with 255 rows should be roughly few MB's.


texts:
I've been importing raster files in ZMAP format using QGIS 2.4 and 2.6. Some of the files work fine but other files for some reason import a black image with nan as a default black and white value. Applying actual min and max values does nothing and shows the values in the raster file as zero.

The ZMAP files have all been exported equally from a gridding package, i.e. same XY min/max, header files look much the same, there are sensible values in all files, etc.

Any suggestions as to what might be causing this?


texts:
I'm trying to catch the OnCreateFeature event throughout the lifespan of the application, but ArcMap seems to be "disabling" my event handler when I exit edit mode in ArcMap. If I subscribe the event handler again after ArcMap "disables" it, it re-enables the first handler and subscribes the second one, leaving with my event handler being invoked twice. How can I keep my events subscribed?

This behaviour is really bizarre, but there's nothing else operating on this event.

Here's how I'm currently subscribing my events.

private void HookupEvents()
{
    var events = (IEditEvents_Event)ArcMap.Editor;
    events.OnCreateFeature += Events_OnChangeFeature;
}

Update

So I found that assigning (IEditEvents_Event)ArcMap.Editor to a static variable stops my events from getting unsubscribed. This leads me to suspect some sort of garbage collection/compiled optimization issue.


texts: How to write code for page setup in Mapbasic

texts:
I am very new to mapbasic scripting. i am creating thematic map using mapbaic code. For that i am searching code for Page setup to set my map as landscape with A4 size. But i didn't found any code for it. does anyone knows it. if any link and resource material is there please reply me.

Thanks and Regards


texts:
I need to merge together 20 32 bit tiff.
For that I tried ArcGis but it keeps crashing (the individual tiff are about 1GB).
I need a finale Tiff in order to run some analysis and I figured that having 20 different file would be messy for the user who is going to do the analysis.

I found a mosaic tool in the raster package in R but I can't really figure out how to use it.
My first idea was to create a list of the tiff and then ask the mosaic tool to process it but it gives me error.

  Error in function (classes, fdef, mtable)  :
    unable to find an inherited method for function "mosaic", for signature "list", "missing"

I found this script online but it seems that many function doesn't work anymore

require(raster)
input.rasters <- lapply(list.files(pattern="^TmB50.*[.]tif$"), raster)
full.extent <- unionExtent(input.rasters)
bounding.raster <- raster(full.extent,crs=projection(input.rasters[[1]]))
res(bounding.raster) <- res(input.rasters[[5]])
resampled.rasters <- lapply(input.rasters, function(input.raster) {
target.raster <- crop(bounding.raster, input.raster)
resample(input.raster, target.raster, method="bilinear")
})
raster.mosaic <- mosaic(resampled.rasters, fun=max)

source http://www.nceas.ucsb.edu/scicomp/usecases/createrasterimagemosaic

Any idea how I  could do that?


texts: Intersecting polygons from different layers based on common values using ArcPy?

texts: I have two polygon shapefiles (ShapeA and ShapeB) that share some areas. The shapefiles have one attribute in common (CODE) and I would like to do a selection, where I select all the polygons in one shapefile one by one and the corresponding polygon in the other shapefile based on same attribute value to get the intersection. Is this possible in some meaningful way using Arcpy?


texts:
I am working on Map Frame layout and would like to remove:
- the oval shape on the edge:  (to make it just 90 degrees)

add  measurment 480 000   490000    and 4520000 450000 corresponding

I tried  Data  Frame Properties -->  not able to fix it.


texts: Passing RasterDataSource object to RasterData object in Javascript API?

texts:
I have a GP service that requires a RasterData object as input. The RasterData object should hold a RasterDataSource object, but I'm not sure of the syntax to connect these two classes. The API reference has no examples. I have tried:

var dataSource = new RasterDataSource();
dataSource.workspaceId = "SharedData";
dataSource.dataSourceName = "raster.tif";
var rasterData = new RasterData();
rasterData = dataSource.toJson();

All properties of the RasterData object return "undefined". What is the correct syntax to do this?


texts: Creating an Atlas where polygons have the same attribute

texts:
I am struggling to create the following Atlas :

I have a coverage layer containing polygons. I want to create an atlas where each map contains polygons that have the same attribute vale (blockname). The image should also be saved as the blockname.

I just cannot figure it out. By the way I am using QGIS 2.6.1. Also I am joining a csv file with attributes to the polygon layer.


texts: Exporting (only) final iteration from ArcGIS ModelBuilder?

texts:
I'm using an iterator in ModelBuilder (1) to merge files, using a feedback loop (2). The output file is named "Output_%n%.txt", (3), so each file is unique.

Now, I'm merging a fair number of files, so I'm trying to speed the process up. Since I don't care about the interim files, I can write them out to in_memory. But I would like to save the final file, viz: "Output_167.txt" to C:\myfile rather than in_memory.

(0) Yes, APPEND would work, but this has other applications.

(1) http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00400000001n000000.htm

(2) http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Iteration_using_feedback

(3) http://resources.arcgis.com/en/help/main/10.1/index.html#//002w0000005s000000


texts: PostGIS SQL error in Java: no spatial operator found for st_contains

texts:
I'm quite new in PostGIS, and I'm trying to run an SQL query. I use this query to obtain a subset of segments from table roads in a group of polygons, then put the result in a separate table.
There are 3 tables:

Table poly: a set of polygons
Table roads: a set of road segments
Table buffer: destination table

The query is like this:
WITH result AS
(
  SELECT r.geom
  FROM roads r JOIN
      (SELECT ST_UNION(geom) area
      FROM poly) p
      ON (ST_CONTAINS(p.area, r.geom))
 )
INSERT INTO buffer
SELECT geom FROM result
;

When I run the query from pgAdmin 4, I can get the result exactly what I want.
However, when I run it from Java, I always get errors:
org.postgresql.util.PSQLException: ERROR: no spatial operator found for 'st_contains': opfamily 16513 type 16400
    at org.postgresql.core.v3.QueryExecutorImpl.receiveErrorResponse(QueryExecutorImpl.java:2553)
    at org.postgresql.core.v3.QueryExecutorImpl.processResults(QueryExecutorImpl.java:2285)
    at org.postgresql.core.v3.QueryExecutorImpl.execute(QueryExecutorImpl.java:323)
    at org.postgresql.jdbc.PgStatement.executeInternal(PgStatement.java:473)
    at org.postgresql.jdbc.PgStatement.execute(PgStatement.java:393)
    at org.postgresql.jdbc.PgStatement.executeWithFlags(PgStatement.java:322)
    at org.postgresql.jdbc.PgStatement.executeCachedSql(PgStatement.java:308)
    at org.postgresql.jdbc.PgStatement.executeWithFlags(PgStatement.java:284)
    at org.postgresql.jdbc.PgStatement.execute(PgStatement.java:279)

I created a Java ANT project in NetBeans 12. I installed the most recent Java, JDBC, Postgresql and PostGIS.
If I execute the common SQL queries without using any spatial function, no problems at all.
Below is my script to run the query.
Class.forName("org.postgresql.Driver");
conn = DriverManager.getConnection("jdbc:postgresql://localhost:5432/gisdb", "user","password");
st = conn.createStatement();
qry= "WITH result AS ( "+
     "SELECT r.geom "+
     "FROM roads r JOIN "+
     "  (SELECT ST_UNION(geom) area FROM poly) p "+
     "ON (ST_CONTAINS(p.area, r.geom))) "+
     "INSERT INTO buffer "+
     "SELECT geom FROM buffer";
st.executeUpdate(qry);

Are there any settings that I missed? I have tried changing JOIN statement, using CTE but no luck. The error remains the same. I never have problems with Oracle or MySQL when calling complex queries with functions from Java.
edit:
I tried to change the username to postgres and the query can run successfully from Java. It seems that my username "user1" cannot access the ST function in public schema. Since I don't want to use administrator username in Java, which privileges should I grant to "user1"?
I tried both grants below with no luck:
GRANT EXECUTE ON ALL FUNCTIONS
IN SCHEMA public
TO user1;

GRANT USAGE
ON SCHEMA PUBLIC
TO user1;


texts: Extracting peaks from DEM and creating 5 m elevation buffer around each peak in ArcGIS Desktop?

texts:
I have a DEM available and would like to identify all the peaks that are likely to be bare rock. I want to be able to isolate these areas and create a 5m elevation buffer around them to represent the non-productive forest. I have used tools like fill but find the output to be to large/general and would like to get more specific..

In summary, I want to find and extract each peak in the DEM to turn it into a polygon out to the closed contour that encompasses the peak. Ideally like the picture attached.


texts: Different min and max NDVI values with different approaches using Google Earth Engine

texts:
I'm trying to calculate the vegetation fraction for land surface temperature calculations of a single Landsat 8 image. I need to find the min and max NDVI of the scene.

I've masked the image and calculated NDVI with two different methods:

var ndvi1 = image.normalizedDifference(['B5','B4']);
var ndvi2 = nir.subtract(red).divide(nir.add(red));

But when I calculate the min and max NDVI values, the results are different for these two methods.

  ndvi1 min: -0.33634447455585664
  ndvi1 max: 1

  ndvi2 min: -94.81818181818181
  ndvi2 max: 63.333333333333336

What is the source of this difference? How should I revise my code for the accurate calculation of min and max NDVI values?

Here is my code:

Code Editor script link

var raw2 = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_174034_20161117');

var vizParams1 = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};
Map.addLayer(raw2, vizParams1, 'raw2');

// Cloud mask of a single Landsat 8 image.
{
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = ee.Number(2).pow(3).int();
  var cloudsBitMask = ee.Number(2).pow(5).int();
  var snowBitMask = ee.Number(2).pow(4).int();
  var waterBitMask = ee.Number(2).pow(2).int();
  // Get the QA band.
  var qa = raw2.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).and(
            qa.bitwiseAnd(cloudsBitMask).eq(0)).and(
            qa.bitwiseAnd(snowBitMask).eq(0)).and(
            qa.bitwiseAnd(waterBitMask).eq(0));
  var image =
      // Scale the data to reflectance and temperature.
      //.select(['B5', 'B7'], ['NIR', 'SWIR']).multiply(0.0001)
      raw2.updateMask(mask);
}

var vizParams2 = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};
Map.addLayer(image, vizParams2, 'Landsat');

print(image,'Landsat');

var nir = image.select(['B5']);
var red = image.select(['B4']);
var ndvi = nir.subtract(red).divide(nir.add(red));
//.rename('NDVI');

var NDVI = image.normalizedDifference(['B5','B4']);

var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};

Map.addLayer(NDVI, ndviParams, 'NDVI image');

Map.addLayer(ndvi, ndviParams, 'ndvi image');

//Calculating min and max values:
var min = ee.Number(NDVI.reduceRegion({
     reducer: ee.Reducer.min(),
     scale: 30,
     maxPixels: 1e9
     }).values().get(0));
  print(min, 'NDVImin');
var max = ee.Number(NDVI.reduceRegion({
      reducer: ee.Reducer.max(),
     scale: 30,
     maxPixels: 1e9
     }).values().get(0));
  print(max, 'NDVImax');

var min = ee.Number(ndvi.reduceRegion({
     reducer: ee.Reducer.min(),
     scale: 30,
     maxPixels: 1e9
     }).values().get(0));
  print(min, 'ndvi-min');
var max = ee.Number(ndvi.reduceRegion({
      reducer: ee.Reducer.max(),
     scale: 30,
     maxPixels: 1e9
     }).values().get(0));
  print(max, 'ndvi-max');


texts: How to fix self intersected polygons visually, as well in QGIS?

texts:
I have same invalid polygon features because of self intersection. I can make them valid by using fix geometries tool in QGIS 3.2. However I want to show them valid visually also. The feature is valid below but visually not, how can I fix it?

It is okay correcting this by using Python and its libraries for instance Fiona, Shapely etc.


texts: Using Tissot's Indicatrix Validate Conversion to ECEF

texts:
When using NED & ECEF as CRS:

https://www.mathworks.com/help/map/ref/ned2ecef.html?requestedDomain=true

Am I should worry from getting 'classical' projection issues?

For example when x,y are too far away from lat0, lon0?

And if so, is Tissot's Indicatrix can help?


texts: How to configure GeoServer WMS with no internet connectivity?

texts:
I am trying to configure GeoServer so that it works with no internet connectivity. When I try to access the GetSchemaExtension, by default, it goes out to the internet for this url:
http://schemas.opengis.net/wms/1.3.0/capabilities_1_3_0.xsd

Is there a configuration out there that allows this file (among many other files) to be configured so that they are accessed internally? I can download the necessary files and upload them onto the server.


texts:
I have the following 9 clusters (each colour represent the appropriate cluster):

I am looking for a clustering algorithm (if any) that would achieve "realistic" results.

The data represents Vancouver's bus stops, crossings, gates, stations etc (as you can see in the left side of the image).

I 've tried K-Means and I obviously had some bad results since it uses the Euclidean distance as its metric and that's not applicable in this scenario. I am using WEKA.


texts: Postgis Self-intersection on geometryraster comparison

texts:
I am rather new to postgis and try to get a (weighted) average height of an area around a point. I am running Postgresql 9.5.5 with Postgis 2.2.1 on Ubuntu. I have imported ASTER GDEM elevation information as a raster and built the following functions, following this tutorial:

WITH statics AS (
    SELECT (ST_DISTANCE(CAST (ST_Project(CAST(ST_SETSRID($1, 4326) AS geography), $2, radians(90.0)) AS geometry), CAST(ST_SETSRID($1, 4326) AS geometry))) AS exp_fact
    -- calculate degrees to expand the point with, $2 is given in meters
)
SELECT SUM(sub.height * sub.weight) / SUM(sub.weight)
FROM (
    SELECT height, ((1 - (dist / $2))) as weight
    FROM get_elev_points_within_geom(ST_MAKEVALID(ST_BUFFER(ST_SETSRID($1, 4326), (SELECT exp_fact from statics))), $1)
    WHERE dist <= $2
) as sub

where get_elev_points_within_geom(geo geometry, center geometry) is:

SELECT ST_Y(ST_CENTROID((a.neighborhood).geom)),
    ST_X(ST_CENTROID((a.neighborhood).geom)),
    (a.neighborhood).val,
    ST_DISTANCE(CAST (ST_CENTROID((a.neighborhood).geom) AS geography), CAST ($2 AS geography)) as dist_m
FROM(
    SELECT ST_INTERSECTION(elevation.rast, $1) AS neighborhood
    FROM elevation
    WHERE ST_INTERSECTS(elevation.rast, $1)
) AS a
ORDER BY dist_m
LIMIT 1000;

Running e.g. this query:

SELECT * FROM get_elev_average_weighted_by_dist_at(ST_MAKEPOINT(8.4305360482003504, 49.013877361340597), 50);

I get:

NOTICE:  Ring Self-intersection at or near point 8.4315277777777791 49.016527777777782
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
SQL function "get_elev_points_within_geom" statement 1
SQL function "get_elev_average_weighted_by_dist_at" statement 1
Total query runtime: 145 msec
1 row retrieved.

I have traced the Problem down to to:

SELECT (ST_INTERSECTION(elevation.rast, ST_SETSRID(ST_GeomFromText('POLYGON((8.4298 49.0131,8.4298 49.0145,8.4312 49.0145,8.4312 49.0131,8.4298 49.0131))'),4326))).geom
FROM elevation
WHERE ST_INTERSECTS(elevation.rast, ST_SETSRID(ST_GeomFromText('POLYGON((8.4298 49.0131,8.4298 49.0145,8.4312 49.0145,8.4312 49.0131,8.4298 49.0131))'),4326))
LIMIT 14;

which will return with the same error message (but different trace, obviously). Interestingly, if I reduce the limit to 13, I will not get the Message.

Also this will run for about 150ms per Query on my machine. This is quite annoying if I query for lots of points.

How can I improve on efficiency?


texts: Wrangling precision in QGIS with a hypothetical planet

texts:
I'm a tabletop roleplaying gamemaster wanting to use QGIS to help organize my notes according to place. GIS addresses my need to store notes and other data that may be anywhere from a few feet from each other to several hundred miles apart, and to store notes that apply to paths and areas. Also, I'm just interested in getting into GIS anyway.
My fictitious planet is a conveniently perfect sphere with a circumference of 25,000 miles. Here's the CRS I cobbled together:
GEODCRS["FictionalPlanet",
    DATUM["Perfect sphere, circumference 25000 miles",
        ELLIPSOID["Perfect sphere, circumference 25000 miles",21008452.5,0,
            LENGTHUNIT["foot",0.3048]]],
    PRIMEM["Prime",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[spherical,3],
        AXIS["geodetic latitude (lat)",north,
            ORDER[1],
            ANGLEUNIT["degree minute second hemisphere",0.0174532925199433]],
        AXIS["geodetic longitude (long)",east,
            ORDER[2],
            ANGLEUNIT["degree minute second hemisphere",0.0174532925199433]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["foot",0.3048]]]

My immediate goal is to set up a grid with an interval of 5ft along the equator and the meridians. With a circumference of 132,000,000 feet that works out to 26,400,000 intervals. In turn, that's 73333.33... intervals per degree. So flip that under 1 and I have 0.00001363636 (last two digits repeating).
If I throw that into the properties window for the grid decoration, I get almost what I want, but not exactly. The distance between two grid lines at the equator and along meridians is 5.027 ft.
I'm guessing this is down to precision limits within QGIS, and ultimately something I will have to live with. But is it possible to improve precision further with some settings somewhere? Or by using a different technique?


texts: How do I fill a Field Combo Box with vector layer fields in QGIS?

texts:
I am new to python and I am trying to write my first plugin for QGIS. Earlier I asked this question to find out how I could populate a combobox with point vector layers (and nothing else).

I've been trying to retrieve the fields of the layers and use them to populate a combobox as well. A couple hours of searching and trying have resulted in the code below. Unfortunately the combobox remains empty even though I am able to get the field names with the code below in the python console in QGIS.

#Fill combobox with vector layer fields.
layerFields = self.dlg.mMapLayerComboBox.currentLayer()

fields = layerFields.pendingFields()
field_names = [field.name() for field in fields]
self.dlg.mFieldComboBox.addItems(field_names)

Code for layer combobox:

layers = QgsMapLayerRegistry.instance().mapLayers().values()
for layer in layers:
 self.dlg.mMapLayerComboBox.addItem(layer.name(), layer)

Answer: my code lacked the setLayer line. Now resulting in the code below:

layerFields = self.dlg.mMapLayerComboBox.currentLayer()
self.dlg.mFieldComboBox.setLayer(layer)
fields = layerFields.pendingFields()
field_names = [field.name() for field in fields]
self.dlg.mFieldComboBox.addItems(field_names)


texts: Find coordinates of extended line segment intersects with given longitude

texts:
I have a line between C1 and C2;

My Coordinates are;

C1: 36.454N - 21.785E
C2: 41.725N - 38.475E

I want to find C3 latitude value based on given longitude when I extend my C1-C2 line segment.

IE C3 Lng is 42.00E.

How can I find it? And I use EPSG:4326 as CRS.


texts:
I have the same problem reported here (QGIS WMS Server), i have one project created in Qgis desktop (tried windows and ubuntu flavors), then exported WMS into Qgis server (ubuntu). All line layers are rendered black loosing their styles. Polygons are rendered ok.

Any ideas?


texts: 2.10 program group option to open QGIS with GRASS

texts:
I have installed QGIS 2.10 on Windows 8.

I cannot find the 'option' to start QGIS with GRASS - as has been previously answered.

Have I possibly missed something in the install? Is this 'missing' from Windows 8


texts: GeoPandas .sjoin polygon in polygon - works fast for some countries, but not for others

texts:
I am new to working with geodata.

Setup: GeoPandas, Python3.6

I have a dataset of square polygons of 10x10km for countries. However, some of these squares are already completely in a neighboring country. To identify and later remove the ones outside the country border.  Hence I have two GeoDataFrames:

squares >> contains all squares plus a country identifier.
countries >>> contains the country borders as polygon or multipolygon

What puzzles me is that the below code works fast for some countries, but takes forever for others. The number of square polygons cannot be the issue as the processing times for countries with approximately 1200 squares can be 45 seconds or 30 minutes.  I suspect that countries with an intricate border (e.g. Sweden and its archipelago) take longer to process.

The code never stops, but it just takes forever. I am thus now running each country one by one rather than the complete list of squares. Each country is saved as a pickle file after processing to clean out memory as it goes up to 7-8GB per country.  Below is my code for Sweden which has some 23,000 squares.

How can I improve speed? I am assigning a spatial index with df.sindex. I also looked at R-Tree as used by Geoff Boeing, but that joins points in polygons. I could not make it work for polygons in polygons. And maybe there is even a better way of doing this.

import pandas as pd
import geopandas as gp
from matplotlib import pyplot as plt
import numpy as np

### Data overview
countries.head(2)
     COUNTRY                                           geometry
0         SE  MULTIPOLYGON (((24.15513 65.81603, 24.12993 65...
1         FI  MULTIPOLYGON (((28.92968 69.05190, 28.82917 69...

# Sample data for Sweden
squares[squares['COUNTRY'] == 'SE'].head(2)
       COUNTRY                                           geometry
717995      SE  POLYGON ((9.64431 58.11107, 9.64343 58.20099, ...
717996      SE  POLYGON ((9.64343 58.20099, 9.64254 58.29091, ...

### Processing
# Subset country
country = 'SE' # SE for Sweden
mask = squares['COUNTRY'] == country
squares_inter = squares.loc[mask, :]
# Spatial join
shape_inter = gp.sjoin(shapeisect, countries, how='inner', op='intersects')
# Amending columns
shape_inter.rename(columns={'COUNTRY_left':'COUNTRY_sq',
                            'COUNTRY_right':'COUNTRY_c'}, inplace=True)
# Get match
shape_inter['match'] = np.where(shape_inter['COUNTRY_sq'] == shape_inter['COUNTRY_c'], 1, np.nan)
# Save to pickle
filename = 'country  + '.pkl'
shape_inter.to_pickle(filename)

### Plotting squares and country borders BEFORE processing
mask = countries['COUNTRY'] == 'SE'
base = countries.loc[mask, ['geometry']].plot(figsize=(10,10), facecolor="none", edgecolor='black')
mask = (squares['COUNTRY'] == 'SE')
squares.loc[mask, ['geometry']].plot(ax=base, facecolor='none', edgecolor='green')

Below is an image that shows how the squares go beyond country borders.


texts: ASP.NET application returns "type “geometry” does not exist" error while Postgis extension is installed

texts:
I have a local Postgres 9.5 with Postgis 2.3.2 installation.
I can confirm Postgis is indeed installed when I watch the extensions, it is installed in the public schema, my application uses a different schema.
The search_path for the database is "$user", public, topology

All the spatial stuff seems to work in the database itself (I have several views that use functions and tables that have spatial data types), however when my application hits something with a spatial type in it, I get the error "type “geometry” does not exist".

The weird thing is that I do not have this problem with my staging and production environment. They run on Postgres 9.4 with Postgis 2.2 installed.
Here the Postgis extension is also allocated to the public schema.

One big difference I found between the two is that 2.3 is not reallocatable anymore. So I cannot set it to a different target schema.

I did a check with a new Postgres 9.4 with 2.2.5 Postgis install (allocated to the public schema at first). I did get the error message, and when I reallocated the postgis extension to the schema my app uses, the error goes away.

How can I fix this for my 9.5 install without reinstalling it entirely?


texts:
I have Navitel Navigator 5.0 map file(.nm3) for navigation.
Is ther is a way to open this file on PC?
Do we have a software that works as a navigator for PC such as Navital Navigation for android.


texts: How to draw a line with OpenLayers by clicking two points?

texts:
I've got website with OpenLayers tile layer. I am trying to draw vector line on the map (only one segment between 2 points) on mouse click and to get coordinates of these two points.  I read a lot on the Internet, but I could not handle it.

That's my code for adding OSM layer, layer switcher, view:

      var view = new ol.View({
      projection: new  ol.proj.Projection({
      code: 'EPSG:3857',
      units: 'm'
      }),
      center: [1030534,5690437],
      zoom: 12
    })

  var layers = [
    new ol.layer.Tile({
      source: new ol.source.OSM(),
      title: 'Open Street Map'
    })
  ];

  var map = new ol.Map({
    layers: layers,
    target: 'map',
    view: view,
    controls: ol.control.defaults().extend([
      new ol.control.OverviewMap()
    ])
  });

var layerSwitcher = new ol.control.LayerSwitcher();
map.addControl(layerSwitcher);

I really have no idea how to make the line-part. I use this
 example, but it's hard for me: https://openlayers.org/en/latest/examples/draw-features.html

I'd like to make something like LineString from the example above but only with one segment.

I tried with that, but it's obviously NOT right. As I said - I'm beginner.

  var source = new ol.source.Vector({wrapX: false});

  var vector = new ol.layer.Vector({
    title: 'Line',
    source: source
  });

     var draw = new ol.Interaction.Draw({
        source: source,
        type: 'LineString'
      });
      map.addInteraction(draw);


texts: How to acquire schema lock is one already exists?

texts:
I'm trying to delete a feature class from a feature dataset or to add a field to it, but I keep getting an error message (pics attached below). I have imported a geometric network I've created from one PGDB to another. However couldn't add any fields. I'm receiving the following error:

Than I planned to delete the feature classes that were imported - I've received the following error:

I have the feeling that both errors are related but couldn't figure out the problem from ESRI help or from other posts in gis.stackexchange.


texts: Geopandas causes recursion error in pyinstaller

texts:
When trying to turn to exe a py file this happens:

RecursionError: maximum recursion depth exceeded

It turns out after a lot of trial and error by deleting libraries that may cause issues that after removing geopandas it finally created the exe.

So the problem now is how to make it work WITH geopandas in it because it is important for my program.


texts: Points in hosted Feature Layer causing rendering issues in ArcGIS Javascript API

texts:
I have Feature Service in my ArcGIS Online account that is made up of a number of feature layers. The feature service was created by publishing a file geodatabase. I am adding the Feature Service as a new GroupLayer() and passing a new PortalItem() to the constructor.

This adds all the feature layers to the map and I can see them in the LayerList. Only features with polyline geometry render properly though. Layers with point geometry cause the GroupLayer to stop rendering (all items in group layer stay at same position on screen even while map is panned/zoomed). It also continuously throws the error Uncaught TypeError: Cannot read property 'getBoundVAO' of null

If I toggle the visibility to false for all of the point layers and leave only the polyline layers then the map behaves exactly as expected.

What is even more bizarre is if I change the visualization properties of the point layers in ArcGIS Online to display as a heatmap instead of points then the rendering issue goes away.

Is there a way to symbolize the point layers as points in the group layer and still have everything render properly? I have a feeling this is a bug in the Javascript API. Can anyone confirm this?


texts: Splitting or decomposing LineString into original single lines or segments using shapely?

texts:
I am trying to split a LineString into single segments:

from shaply import LineString,MultiPoint
from shapely.ops import split
xp = [-2,1,2,2.5,4,1,-1,2,4,5,6]
yp = [-2,1,1,0,1,2,2,3,4,4,4]
list_of_point = list(zip(xp,yp))
profile = LineString(list_of_point)

At the moment I am using split:

points_to_split = MultiPoint([Point(x,y) for x,y in profile.coords[1:]])
splitted = split(profile,points_to_split)
for i in splitted:print(i)

I've got the expected result but it seems quite overcomplicating a process that should be quite easy to handle internally.

Am I missing something?


texts: Collector - Related tables and default filling of field

texts:
Have a polygon feature class with a related table in collector.
Is there a way for when a new record is added in the table for a field (area) to be automatically calculated from the features geometry?
The features will be edited over time and I want to be able to store the current area in the table of when the record was created to perform a calculation later on.


texts: What is the difference between "due North", "N0degreesW" and "N0degreesE"?

texts:
Is there a difference between "N0°W" and "N0°E"?

These terms have been used in connection with Strike Lines.

The final bearing computed was exactly 0.

Will the direction East and West matter if one is chosen over the other or if "due North" or simply "North" is chosen instead?


texts: Update record based on series of conditionals using Python

texts:
I have a feature class which I've been using a relate to compare duplicate field values and was thinking of a way to eliminate the manual process of checking these duplicate values and making a decision as to which records to keep. There are several parameters which determine which record is kept vs. which one is discarded. For now I'd like to flag which records are kept so I can do some spot checking on the attributes.

I have some mock data below to show the desired results in a "Flagged" field. The conditions would involve first using a where clause to prompt the cursor to scan the "ID" field where ID = X. For this selection compare fieldB records, if they are equal values and one of them does not contain the letter 'C'
compare fieldA records and flag the one with the higheest numeric value unless
there is a null in one of the fieldA records. If one of these records contains a null value, compare the records in fieldC and flag the greatest value.

Feature class records:

fieldA  fieldB  fieldC  Flagged ID
1       ABC     4                X
2       AB      5                X
3       A       6                X
1       EFG     4                Y
2       EF      5                Y
3       E       6                Y

Desired result in "Flagged" field:

fieldA  fieldb  fieldC  Flagged ID
1       ABC     4                X
2       AB      5                X
3       A       6         1X     X
1       EFG     4                Y
2       EF      5                Y
3       E       6         1Y     Y


texts: Create kmz with coordinates from images and thumbnails for symbols

texts:
I have about 110 geotagged images and am looking for an effective way to create a .kmz which contains for each image a point based on its geotag and uses a thumbnail of the picture as symbol. Clicking on the thumbnail should - as usual - open the larger image.

I am able to manually load all images into Google Earth, but this takes a long time and I am not able to create the thumbnails.


texts: Processing daily PRISM climate Raster datasets?

texts:
I am looking to perform principal component analysis (PCA), based on daily gridded raster datasets (7 different climate variables) between 1985 and 2015.  I was hoping to take each daily dataset (a single raster) and create an average for, say, an entire year, or even for a season of time.  However, it would take a lot of work and time, using the Raster Calculator tool, to come up with an average raster for that length of time.

Would it be possible to automate a task like this, using ModelBuilder within ArcMap?


texts: QGIS 2.12.3 Crashes on intersecting two vector layers

texts: I'm trying to intersect two vector layers (same projection, no spaces or unallowed characters in the filename) using intersect tool, but qgis crashes at every try. It says something like "created minidump in C:" and shows a path to a dmp-file. I already uninstalled Pisa (2.10 Version) and installed the Lyon-Version, but I still have the same problem. Can anybody give some hint to solve this problem?


texts: Unable to process multiple polygons using SpatialPolygons object

texts:
I have a list with 4 objects, each object is a data frame. The data frame is the coordinate matrix of each polygon. I would like calculate the total area of the 4 polygons and then draw the voronoi tessellation.
Here is my workflow,
first, I convert the data frame within the list to polygon:
  Region_HE_srl <- lapply(Region_HE,FUN= function(x) Polygon(x))

Then, I create the SpatialPolygons object:
  Region_HE_sp <- SpatialPolygons(list(Polygons(Region_HE_srl,length(Region_HE))))

Now that the SpatialPolygons object is created:

I used the 'gArea' function from the 'rgeos' package to calculate the area.
  area <- gArea(Region_HE_sp, byid = TRUE)

Normally, when there is only 1 polygon within the SpatialPolygon object, this function works fine, but now as there are 4 polygons, an error occurs:

Error in createPolygonsComment(p) :
rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 1

Same thing when I applied the voronoi.polygons function from 'SDraw' package.
  voronoi.tess <- SDraw::voronoi.polygons(pts_ppp, Region_HE_sp)

The error occurs again. As the package claimed that the bounding polygon parameter should be a spatialpolygons object, I cannot figure out why this error happened.
Any suggestions?


texts: How to calculate average NDVI values for buffers in ArcGIS?

texts: I would need to calculate an average NDVI value for each buffer among everal hundred buffer-areas. The buffer areas are built around observation points. I have all buffer areas in one .shp file and NDVI values in a 32 Bit floating point .tif file. I have also the Band 3 and Band 4 values in separate .tif files (if that helps). I am using ArcMap 10.2 and have e.g. the Spatial Analyst Tools package. What should I do to get the average NDVI values for each buffer?


texts: Why do multiple raster processing tools fail with this DEM?

texts:
I have a provincial sized DEM. I need to clip it to a fairly large watershed boundary. Should be easy. I've tried:

Data Management Tools -> Raster -> Raster Processing -> Clip

(Did I make sure "Use Input features for clipping geometry and "maintain clipping extent" where checked? Yup). Here's the error messages:

Says "No spatial reference exists" but it's a liar:

Extract by Mask under Spatial Analyst tools doesn't work either.

I thought maybe my boundary shapefile was the problem. Created a new one using the draw tool, makes no difference. Shapefile set to same datum and projection as DEM.

The DEM itself won't project into a different coordinate system. It won't export - I've tried exporting to tiff, the png, to grid... an error always occurs. It is an unworkable DEM and I do not know why. I've clipped DEMs before no problem.

Can you see anything about this DEM that might make it behave this way? Here's it's properties.


texts:
I would precess MODIS images on QGIS (in particular swath images MOD09A1).
I use MODIS Reprojected Tool Swath to trasform .hdf file in GEOtiff but there are different avaible bands (16) and i don't know what are the reflectance bands number 2 and 7. I would process these band in QGIS to elaborate NDII index.


texts: Importing area attributes from one shapefile into another grid shapefile?

texts:
I have two shapefiles of same extent. One is a shapefile of tree canopy coverage with the area attribute. The other shapefile is 20m*20m grids of the same area, now I have to import the canopy coverage information into the grid shapefile attribute. Note that the canopy coverage is falling in multiple grids but i need to clip out the information falling exactly within the grids.

I have tried with spatial join but the issue is with the canopy sizes the tool uses intersection which is selecting whole canopy layer which is falling in multiple grids. I need to get the clipped area information within each grid.


texts:
I've got a postgresql table that has the following attributes:

name        String
begintime   String
endtime     String
the_geom    Point

Now I'd like to create an sql view in geoserver (2.2 snapshot). In short, i only want to map the features where begintime equals a certain param. Looking at the geoserver documentation and the questions here at the forum, I came up with this:

select name, begintime, endtime, the_geom
from n80_tellingen_merge
where begintime = %utct%

This gets me the following error:

ERROR: syntax error at or near ")" Position: 111

So basically I'm asking two things:

What am I doing wrong?
Since I don't have brackets or more than 111 positions, where can I find the actual query sent to postgresql? Can I consult a log somewhere? (Ubuntu server 13.04)

My apolagies if I'm missing something obvious.


texts:
I'm using QGIS 2.14 and GDAL 1.11.3 and I'm running into a weird error with the clipped outputs.  I created a shapefile in QGIS and used it to clip my image.  Here's a simplified version of the command the Clipper plugin generates:

gdalwarp -q
         -cutline ${TESTDIR}/shapefile.shp
         -crop_to_cutline
         -tr 5.0 5.0
         -of GTiff
         ${TESTDIR}/${INFILE}.TIF
         ${TESTDIR}/${INFILE}_sub.tif

My results are below. The shapefile is displayed in the dotted red line, with the clip overlaid on the original.  Notice the grid offset produced by gdalwarp.  I've also run this on lower resolultion Landsat8 data with similar results (the ratio of the offset error is the same, i.e. ~20% in x, ~5% in y).
This issue QGIS clipper problem when mask used is similar except I am already using the -tr option.  Anyone know what's going on?  Is this a GDAL bug or a QGIS?


texts:
I have a spatialite database I created that spans several years of crime data. To facilitate a real-world scenario where I may simplify the interaction for a user by introducing a View that simply limits the time span or what types of crimes are returned, I want to know if it is possible for QGIS to read this View like it would any other spatial table in my database. The point here is that I obviously do not want to create a separate table for every possible View I may use, and I am trying to make it simple by not requiring the user to know the database and the appropriate SQL to subset the data themselves. A similar workflow in the ArcGIS world might be creating separate layer files that the user could open in ArcMap that access the database (say, SQL Server) and present only the limited data set.

I created a View in my database, but it does not show up in the list of feature classes I can import when I connect to the spatialite db in QGIS. Maybe Views aren't considered spatial? (This is something I'm going to look into later, but I thought I already did; I just don't recall the conclusion.) In any case, how might one replicate the sort of work flow I want or the similar ArcGIS alternative?


texts:
How do I specify the polygons inside a circle.

I mean the polygons inside the circle completely and also the polygons that partly are inside the circles.


texts:
I want to create SLD with an external graphic in GeoServer, I used this:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http**strong text**://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>City Point</Name>
    <UserStyle>
      <Title>City Simple Points</Title>
      <FeatureTypeStyle>
        <Rule>
        <Title>city</Title>
            <PointSymbolizer>
             <Graphic>
              <ExternalGraphic>
                <OnlineResource xlink:type="simple" xlink:href="city.svg"/>
                <Format>image/svg+xml</Format>
              </ExternalGraphic>
              <Size>
                <ogc:Literal>20</ogc:Literal>
              </Size>
             </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

I saved the svg file in this path: /var/lib/tomcat8/webapps/geoserver/data/styles/ but it doesn't work.


texts: Variable substitution in SLD for geoserver 2.13 RasterSymbolizer ChannelSelection

texts:
I'm trying to publish a layer with a style that allows the extended WMS request to specify what image bands to map to the RGB values. Based on existing examples I assumed that my SLD should have something like the following RasterSymbolizer element:

<RasterSymbolizer>
        <ChannelSelection>
          <RedChannel>
            <SourceChannelName>
              <ogc:Function name="env">
               <ogc:Literal>redchan</ogc:Literal>
               <ogc:Literal>1</ogc:Literal>
              </ogc:Function>
            </SourceChannelName>
          </RedChannel>...

But the SLD validation fails with the following error:

Element 'SourceChannelName' is a simple type, so it must have no element information item [children]

I am following this documentation:

http://docs.geoserver.org/stable/en/user/styling/sld/reference/rastersymbolizer.html
http://docs.geoserver.org/stable/en/user/styling/sld/extensions/substitution.html

Is there a way to do this or is it not supported for the ChannelSelection element? Are there any more examples and/or documentation for SLD variable substitution that I can look at?


texts:
I have a DEM of a riverbed and i want to know if the elevation can be increased by different values on the left and right side of it and in between them the values to increase with an average of the two values as shown in the photo.
Basically, on the left side, the height will be increased by 25 from 10 to 35 while on the right the values will be increased by 35 from 10 to 45. The values in between them will be increased with a linear average so basically, the middle value of 15 will be increased by 30 and will get to 45. Is there any way to do this?


texts: What does 'projected bounds' and 'WGS84 bounds' mean in the definition of a CRS (on spatialreference.org)?

texts:
The site http://spatialreference.org gives 'Projected Bounds' and 'WGS84 Bounds' for each Coordinate System. What does these 'bounds' mean? e.g. for ETRS89 / UTM zone 32 N this is given:

I know what the difference between 'WGS84' and 'Projected' means. I want to know why there are 'bounds'.

In epgs.org there are no such bounds. For 'ETRS89 / UTM zone 32 N' epsg.org says:

So epsg.org includes Germany (which is valid) while the bounds of spatialreference.org did not include Germany.

(I want to understand what the 'bounds' mean because AutoCAD Map 3D seams do use these bound when consuming a WGS service.)