texts: ESRI UK Productivity Suite: Critical Error - Help Resolving

texts:
I am using the OS Conversion tool under ProductivitySuite v2.1. to convert OS MasterMap gml to a feature class.

Having carefully followed the instructions I am getting the following: Error 18/05/2012 09:56:57 CRITICAL ERROR: Unable to validate OrderNumber

This is the first time I have used the tool and unsure as to How would I resolve this?

Apparently the data has been copied directly from the CD as it was provided by the OS.

Log File below, if it helps:

18/05/2012 09:56:57 Log File created.
Info 18/05/2012 09:56:57 =================================================
Info 18/05/2012 09:56:57 OS Data Converter
Info 18/05/2012 09:56:57 Version: 2.1 Build: 557
Info 18/05/2012 09:56:57
Info 18/05/2012 09:56:57 Start Processing
Info 18/05/2012 09:56:57  Profile name: OS MM SWW
Info 18/05/2012 09:56:57  Start condition: New
Info 18/05/2012 09:56:57
Info 18/05/2012 09:56:57 =================================================
Info 18/05/2012 09:56:57 Validation:
Info 18/05/2012 09:56:57  - Validating data order number
Info 18/05/2012 09:56:57
Error 18/05/2012 09:56:57 CRITICAL ERROR: Unable to validate OrderNumber
Error 18/05/2012 09:56:57 CRITICAL ERROR DETAILS:
Error 18/05/2012 09:56:57   ESRIUK.Products.UKPS.DC.Core.DataConversionException: Unable to validate OrderNumber ---> System.Exception: Unable to validate order number
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Profile.OrderNumber.Validate(List1 sourceItems)
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Profile.OrderNumber..ctor(List1 sourceItems)
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Profile.DataConverterSourceList.ValidateOrderNumber()
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Core.DC.DataConverterProcess.ValidateOrderNumber()
Error 18/05/2012 09:56:57      --- End of inner exception stack trace ---
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Core.DC.DataConverterProcess.ValidateOrderNumber()
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Core.DC.DataConverterProcess.DoValidate(String paramName)
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Core.DataConverterEngine.ValidateProfile()
Error 18/05/2012 09:56:57      at ESRIUK.Products.UKPS.DC.Core.BaseEngine`3.ProcessProfile(T_Profile profile)
Info 18/05/2012 09:56:57
Info 18/05/2012 09:56:57 OS Data Converter

many thanks....


texts: Spatial indexing of a point pattern in R

texts:
I have a point pattern and a tessellation (grid) in R. Is there a way to spatially index each of my points? In other words, my tessellation grid has n quadrats, I would like to know for each of my points which of these (1:n) quadrats they are associated with.
If not, is there another way to do so in R?
Many thanks for your suggestions!


texts:
I'm on QGIS 3 and have used the 'AusMap' plugin to upload google maps raster images and a 5m LiDAR DEM. I have created a shapefile (polygon) and now I want to clip the raster files to the extent of the polygon. Whether I 'clip by extent' or 'clip by mask layer', I still run into the same issue--an error 4.
Here is the error message:

"2021-02-24T22:56:39     INFO    gdal_translate -projwin 151.917322934
-32.444796101 151.936019747 -32.45118951 -of GTiff "identifier=1&url=http://services.ga.gov.au/gis/services/DEM_LiDAR_5m/MapServer/WCSServer?"
"C:/Users/curti/OneDrive/Documents/GIS Base/the_block/clip.tif"
2021-02-24T22:56:39     INFO    GDAL execution console output
ERROR 4: `identifier=1&url=http://services.ga.gov.au/gis/services/DEM_LiDAR_5m/MapServer/WCSServer?'
does not exist in the file system, and is not recognized as a
supported dataset name."

There seems to be some issue with the LiDAR file, but I don't understand what the problem is or how to solve it.


texts: Seeking Chesapeake Bay Watershed Boundary Shapefile?

texts:
Does anyone have the shapefile that depicts the boundary of the Chesapeake Bay Watershed?

I have looked through countless GDB's from government agencies without any success.


texts: Saving Result File from command line ArcPy?

texts:
What do I need to do in order to get a custom arcpy script to save an execution run into the Results window, and/or just straight to an .rlt file directly?

The script is run from the command line and isn't in a toolbox.

Motivation: I use results as documentation, to see what has been done in a particular project, and to recreate particular steps dynamically at will, confident that I will get the same outcome.  Significant steps are saved as .rlt for future reference and re-use in other projects. I want to capture my non-toolbox scripts as well as model scripts in this manner as well as stock tools.


texts: GeoTIFF is stretched sideway when translating CRS to EPSG:4326 from EPSG:32654

texts:
I want to convert a GeoTIFF's coordinate reference system to EPSG:4326 from EPSG:32654. But the output image is stretched to the side.
I have tried two different methods (gdalwarp and rasterio), but the result is the same.
Why does this happen? If you know of a solution, I would like to know how to solve this problem.
This is the command I used.
gdalwarp -s_srs epsg:32654 -t_srs epsg:4326 source.tif target.tif

Here is my GeoTIFF file EPSG:32654 (before conversion) and its gdalinfo.

Driver: GTiff/GeoTIFF
Files: source.tif
Size is 9105, 5977
Coordinate System is:
PROJCS["WGS 84 / UTM zone 54N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",141],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32654"]]
Origin = (286156.478620000008959,3945902.046580000314862)
Pixel Size = (0.018810000000000,-0.018810000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  286156.479, 3945902.047) (138d38'18.28"E, 35d38' 1.71"N)
Lower Left  (  286156.479, 3945789.619) (138d38'18.39"E, 35d37'58.06"N)
Upper Right (  286327.744, 3945902.047) (138d38'25.09"E, 35d38' 1.84"N)
Lower Right (  286327.744, 3945789.619) (138d38'25.20"E, 35d37'58.20"N)
Center      (  286242.111, 3945845.833) (138d38'21.74"E, 35d37'59.95"N)
Band 1 Block=9105x1 Type=Byte, ColorInterp=Red
  NoData Value=-10000
Band 2 Block=9105x1 Type=Byte, ColorInterp=Green
  NoData Value=-10000
Band 3 Block=9105x1 Type=Byte, ColorInterp=Blue
  NoData Value=-10000
Band 4 Block=9105x1 Type=Byte, ColorInterp=Alpha
  NoData Value=-10000

Here is my GeoTiff file epsg:4326 (after conversion) and its gdalinfo.

Driver: GTiff/GeoTIFF
Files: target.tif
Size is 9709, 5310
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (138.638412488183292,35.633844968775641)
Pixel Size = (0.000000197720465,-0.000000197720465)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 138.6384125,  35.6338450) (138d38'18.28"E, 35d38' 1.84"N)
Lower Left  ( 138.6384125,  35.6327951) (138d38'18.28"E, 35d37'58.06"N)
Upper Right ( 138.6403322,  35.6338450) (138d38'25.20"E, 35d38' 1.84"N)
Lower Right ( 138.6403322,  35.6327951) (138d38'25.20"E, 35d37'58.06"N)
Center      ( 138.6393723,  35.6333200) (138d38'21.74"E, 35d37'59.95"N)
Band 1 Block=9709x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=9709x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=9709x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=9709x1 Type=Byte, ColorInterp=Alpha


texts: QGIS Voronoi polygons being generated incorrectly

texts:
I am trying to generate some Voronoi/Theissen polygons from points across the UK. From my understanding, any point within a Voronoi polygon should be closer to that polygons origin point than any other origin point; however, this isn't whats happening. The most obvious example is below:

In this example the Swansea (the bay which is close to the centre point of my measuring line) is about 33km closer to the Pembrokeshire point than the Devon point, but that region is still covered by the Devon Voronoi polygon.
Am I setting this task up incorrectly? I have since tried to use GRASS's function which yields the same result, SAGAs which I can't get working, and MMQGIS's function, which also yields the same result. Or is this working correctly, and I'm misunderstanding the purpose of these tools/Voronoi diagrams?


texts: How to get the regional maximum or minimum of raster data by using PostGIS?

texts: If I have Mean sea surface pressure data in raster format, is it possible for PostGIS to look for regional maximum or minimum? So that I can locate the L or H as the chart shown below? (image source: NOAA)


texts: Increased memory usage while zooming and dragging

texts:
I have a problem that deals with memory. When I run my GIS program, it uses 495K from my memory after 1 minute. And it continues to work with same value.

But when I zoom in / zoom out / drag the map to other sides, the memory usage of my programme increases to 500-600K. and it does not decrease back to old value, it becomes constant value after the zooming and dragging stops.

I am clearing the unused memory all 1 second with:

GC.AddMemoryPressure(GC.GetTotalMemory(true));

GC.Collect(GC.MaxGeneration, GCCollectionMode.Forced);

GC.Collect();

How can I resolve this problem?


texts: Convert a geotiff web-mercator into an azimuthal equidistant projection

texts:
I have a geo-referenced tiff file, in the web-mercator projection; and I need it converted to an azimuthal equidistant projection.
I'm new to GIS and can't figure out how to do this.

I'm open to using any program for any OS, as long as it gets the job done.

And as far as OS-specific is concerned, I have computers with : Linux-Debian, Mac OX 10.6.8, Windows 7, and Windows XP.


texts:
I have .shp file of polygons and single (point) location.

How can I retrieve the closest polygon to that single location?


texts: What limits Google Maps from displaying when zooming in tighter than 1:1,411 in QGIS?

texts:
In QGIS 2.18.15 while using Google Maps as a basemap anytime I zoom in tighter than 1:1,411 Google Maps stops displaying.

What is limiting Google Maps from displaying when attempting to zoom in tighter?

I've been searching for an answer, but I don't seem to find anything.  I've installed QGIS as 32bit and I don't have (never did) the "Tile Map Scale Plugin" installed.


texts: QGIS2web page does not display parallel and curved labels to lines

texts: I'm working on a web map with qgis 2 web. But I have a problem. When I export the map the street names do not appear in parallel with the lines but horizontally.


texts:
I want to upload a jpeg or some raster image into a C# form and then on entering the latitude and longitudes I need to be able to create an GEOTIFF file of the same image. How can this be done in C#.NET?

you had asked me to use GDAL libraries i did add references,
i'm using Visual studio 2010, i have included library files,
i have computed llx lly lux luy also, but how to tag it to an image, is there any inbuilt function or anything?

this is the code i have used

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using OSGeo.GDAL;
using OSGeo.OGR;
using OSGeo.OSR;

namespace geo_tagging
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }

        private void button2_Click(object sender, EventArgs e)
        {
            int orx = pictureBox1.Image.Height;
            int ory = pictureBox1.Image.Height;
            int ctrx, ctry, lux, luy, llx, lly, bufx, bufy;
            bufx = orx / 108000;
            bufy = ory / 108000;
            ctrx = int.Parse(textBox2.Text);
            ctry = int.Parse(textBox3.Text);
            lux = ctrx - bufx;
            luy = ctry - bufy;
            llx = ctrx + bufx;
            lly = ctry + bufy;

            Gdal.GCPsToGeoTransform(
        }

        private void button1_Click(object sender, EventArgs e)

        {
            OpenFileDialog dlg = new OpenFileDialog();

            if (dlg.ShowDialog() == DialogResult.OK)
            {
            PictureBox PictureBox1 = new PictureBox();
            string Chosen_File = "";
            Chosen_File = dlg.FileName;
            textBox1.Text = dlg.FileName;
            pictureBox1.Image = Image.FromFile(Chosen_File);
            int X = pictureBox1.Image.Height;
            int Y = pictureBox1.Image.Width;
            label2.Text = X.ToString();
            label4.Text = Y.ToString();

            }

            dlg.Dispose();
        }

        private void button3_Click(object sender, EventArgs e)
        {
            textBox1.Text="";
            textBox2.Text="";
            textBox3.Text="";
            pictureBox1.Image.Dispose();

        }

       private void pictureBox1_Click(object sender, EventArgs e)
        {
            OpenFileDialog dlg = new OpenFileDialog();

            if (dlg.ShowDialog() == DialogResult.OK)
            {
                PictureBox PictureBox1 = new PictureBox();
                string Chosen_File = "";
                Chosen_File = dlg.FileName;
                textBox1.Text = dlg.FileName;
                pictureBox1.Image = Image.FromFile(Chosen_File);
                int X = pictureBox1.Image.Height;
                int Y = pictureBox1.Image.Width;
                label2.Text = X.ToString();
                label4.Text = Y.ToString();

            }

            dlg.Dispose();
        }

       private void textBox1_Click(object sender, EventArgs e)
       {
           textBox1.Text = "";
       }

        }
 }

never mind got it guys, we need to use a exe file to call it as a process and pass arguments into it!

thanks for the help


texts:
I have a shapefile (.shp) of 180k roads with information of street name.

I need to do a Dissolve on that road network to make only one feature per name.

But I have a problem by doing that task because I have cases that the name differs only by one character.

For example:

Rua do servidor pÚblico & Rua do servidor pUblico

or

Rua do SouZa & Rua do SouSa

I was thinking to do something like select all the attributes that match 99% of the string to find that cases.

I don't know how to do it, maybe a python or VB script will do it.
Please, give me a light, and I'll follow it.

I've found something like FUZZY MATCHING...I'll research about it.


texts:
I would like to create a virtual raster from an .img file that has associated .ige, .rde, .rrd, and .xml files. I am doing this to enable virtual I/O so I want an "exact copy" of the original raster.
I am just calling this from the command line:
gdalbuildvrt file.vrt file.img
This creates a VRT but the attribute table is stripped; instead of categorical data I now have continuous values. How can I fix this?


texts: ArcObjects - Adding ArcGIS Server Layer to MXD Programtically

texts:
I am trying to add a layer from ArcGIS online to a mxd programatically with no luck.  See code sample below.  The error is happening at pAGSServerConnectionFactory.Open(props, 0);

Anyone been able to do this?

private void AddArcGISOnlineLayer(IMap map)
        {
            IAGSServerConnectionName servConnName = new AGSServerConnectionNameClass();
            IPropertySet props = new PropertySetClass();
            props.SetProperty("URL", "http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer");
            servConnName.ConnectionProperties = props;
            IAGSServerConnectionFactory pAGSServerConnectionFactory = new AGSServerConnectionFactoryClass();
            IAGSServerConnection pAGSConnection = pAGSServerConnectionFactory.Open(props, 0);
            IAGSEnumServerObjectName pEnumSOName = pAGSConnection.ServerObjectNames;
            IAGSServerObjectName pSOName = pEnumSOName.Next();

            while (pSOName != null)
            {
                if (pSOName.Name == "Layers")
                    break;
                pSOName = pEnumSOName.Next();
            }

            IName pName = (IName)pSOName;
            IMapServer mapServer = (IMapServer)pName.Open();
            IMapServerLayer msLyr = new MapServerLayerClass();
            msLyr.ServerConnect(pSOName, mapServer.DefaultMapName);
            IMapServerGroupLayer group = (IMapServerGroupLayer)msLyr;
            ILayer msLayer = (ILayer)msLyr;
            map.AddLayer(msLayer);

        }


texts: Geopandas tick labels from a diferent coordinate system

texts:
I have a GeoDataFrame that includes points with the relevant geometries (coordinates system WGS84 EPSG:4326)

gdf:
Name      x     y   Coordinates
E   -4.375  51.5    POINT (-4.375000000000008 51.49999999999999)
N   -5.000  52.0    POINT (-4.999999999999994 52)
NE  -4.375  52.0    POINT (-4.375000000000008 52)
NW  -5.625  52.0    POINT (-5.625000000000006 52)
S   -5.000  51.0    POINT (-4.999999999999994 50.99999999999999)
SE  -4.375  51.0    POINT (-4.375000000000008 50.99999999999999)
SW  -5.625  51.0    POINT (-5.625000000000006 50.99999999999999)
W   -5.625  51.5    POINT (-5.625000000000006 51.49999999999999)

I would like to plot with a background map tiles served from openstreetmap. I understand in order to do that in need to change the coordinate system to Pseudo-Mercator WGS84 EPSG:3857. Then request the tile with contextily.

gdf = gdf.to_crs(epsg=3857)
xmin, ymin, xmax, ymax = gdf.total_bounds
url = 'http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'
basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=8, url=url)

And then plot the basmap and the points

ax = gdf.plot(figsize=(10, 10))
ax.imshow(basemap, extent=extent, interpolation='bilinear')

This works fine but I would like the ticks labels to be of the original coordinate system (EPSG:4326). I was thinking that I could get_xticklabels and get_yticklabels and then do coordinate transformation with gdal but there mast be an easier way.


texts: How can I integrate railway embankments to a heightmap?

texts:
I would like to create a 16bit PNG heightmap for a PC game, which simulates a real landscape, and I have created a QGIS project.
In this project there is a raster layer for the terrain, which is based on an "*.hgt"-file (30m resolution).
There is also a vector layer that shows the railroad tracks over the landscape. I added this vector layer to the QGIS project using the Quick OSM plug-in (key: railway, value: rail).
I would now like to merge both layers in a heightmap / in a GeoTIFF in such a way that the tracks are integrated on the terrain as slight elevations or depressions - in the same way that they run over railway embankments or in terrain incisions in the real world.
This works probably only if the tracks in the vector layer contain height information - if not, I would like to integrate all tracks as slight increases (1-5 meters) on the terrain, so that railway embankments are automatically generated as part of the heightmap.
How can I implement this?


texts: ArcPy CAD to File Geodatabase conversion error 000732?

texts:
I have 19 unprojected CAD files in a Folder. Some are .dwg and .dxf extension. I am attempting to convert these CAD files to a GDB using the arcpy.CADToGeodatabase_conversion() tool

I wrote this little script below and it gives me an error

import arcpy
import os
CADpath = 'path\\CADstuff'
GDB = 'path\\CAD.gdb'
sr = '3424'
for cad in os.listdir(CADpath):
    cadname = str(cad[:-4])
    arcpy.CADToGeodatabase_conversion(cad, GDB, cadname,'',sr)
    print cadname

  Traceback (most recent call last):
  File "path\CADtoGDB.py", line 9, in <module>
    arcpy.CADToGeodatabase_conversion(cad, GDB, cadname,reference_scale,sr)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\conversion.py", line 1637, in CADToGeodatabase
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input CAD Datasets: Dataset 080401_DLU_2012-06-21.dwg does not exist or is not supported
Failed to execute (CADToGeodatabase).

*GDB created before script
*First time working with CAD files


texts:
I've been experimenting with ArcGIS Pro's "Visualizing Space Time Cube" tool but found that there is a feature limitation when I try to publish it via ArcGIS Online. My goal is to publish it on the web so end-users can interact with the cube for the whole study area. However, using ArcGIS Pro limits the study area which leads me to exploring the open source way.

I did some searching but my initial findings came up short. Has there been a FOSS solution in creating a space time cube? The subject matter is mapping fire occurrences of a national park over 15 years using a bin size of 3 months/1 mile.


texts: How to extract values based on names of rasters?

texts:
I have a raster r that contains values for each grid cell (which vary from one cell to another). I have other 24 rasters where the names of the rasters are:

the name of r1 is data-00.img,
the name of r2 is data-01.img,
.........
the name of r24 is data-23.img

What I want to do is to look at the first pixel in r if NA,return NA but if there is value for example 12, search in the names of the other rasters (should be data-12.img) and then extract the corresponding value and put it instead of the 12 value in r.finally, do the same for all pixels.

 r <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
 r1 <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
 r2 <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))

My real data r could read as:

      cod <- file("C:\\Users\\data.img","rb")
      ne<- readBin(cod, numeric(), size=4,  n=1300*500, signed=TRUE)
      r <-raster(t(matrix((data= cet), ncol=500, nrow=1300)))

the other 24 rasters could be read as:

       dir2<- list.files("C:\\Users\\datasm", "*.img", full.names = TRUE)
       st <- stack(dir2)


texts: Cutting Census Block Data to a Study Area

texts: I'm trying to cut census block data to a study area. I have a study area polygon (made from selecting whole counties of interest) and clipped the census data by this. When I visited the new Census file I found this process has created "slivers" from the original census dataset that are very small (area is less than a mile) but its population is in the thousands; the entire population of the original block has been cramped into this tiny area. I had thought clipping the data by the study are would've been successful since block groups, at least to my knowledge, do not cross counties. Is there any way to easily remove these extraneous slivers or a better way to cut the census data to a study area?


texts: How do I compare two features to test if they are identical?

texts:
What is the most efficient method for determining if two features (points, lines, or polygons) are geometrically identical? I'm not concerned about a square figure that contains eight vertices (four corners and vertices at segment midpoints) being flagged as identical to a "true" square (four corners and no vertices at segment midpoints).

Potential gotchas are the order of the vertices and direction of line or polygon. Suppose one square is defined starting in the northwest corner and another square starts in the southeast. Coordinates of the vertices and all other properties are identical. I would expect these to be a match since just the order of vertices is just different. Or what if one square is defined clockwise but the other counter-clockwise? Again, this shouldn't matter.

What is the most efficient way to find if the two are identical?


texts: Looking for the "inverse" of the NEAR tool

texts:
To calculate wind fetch over water I need to create a point on the shoreline and then calculate the distance to to other shorelines all directions (0 - 360 deg). The shortest distance to another land feature = the short fetch = low waves. On the other hand, a long fetch or distance to another shoreline feature allows large waves to propagate. (In the extreme case fetch may be huge- crossing the Gulf of Mexico for example- in which case fetch is considered infinite).

How can I perform this calculation? The NEAR tool does a great job at measuring the distance to the nearest feature but not the furthest. This seems as though it should be simple but I haven't been able to come up with a clean workflow.


texts: Is it possible to query raster from ArcGIS Server 10.1 REST?

texts: I have a map service that contains a SDE raster layer and served out by ArcGIS Server 10.1. Is it possible and how to query through REST service on the raster layer? Basically what I want to do is to query the pixel value for the given x, y coordinates. Didn't find any example to do it.


texts: How to select island (donut hole) polygons witin the same polygon layer?

texts: i have a dataset containing many polygons. In this dataset some polygons have other polygons (as a donut) inside of them. So those island polgons are contained by only one polygon but it originated from the same layer and has the same attributes. I want to select and merge them but i dont know how to select/find them. It is possible to select polygons based on the number of donuts inside of them?


texts: Convert NIFC GeoJSON daily data to KML file using GeoTools

texts:
I am not able to convert the daily NIFC GeoJSON data into KML file using GeoTools lib. I am pretty new to GIS data.
Sample GeoJSON data : https://opendata.arcgis.com/datasets/5da472c6d27b4b67970acc7b5044c862_0.geojson
I have tried the below code so far:
URL url = URLs.fileToUrl(new
 File("C:\\datasets\\5da472c6d27b4b67970acc7b5044c862_0.geojson"));
HashMap<String, Object> params = new HashMap<>();
params.put("url", url);
DataStore in = DataStoreFinder.getDataStore(params);

if(in == null) {
  throw new IOException("couldn't open datastore from "+url);
}
SimpleFeatureCollection features = in.getFeatureSource(in.getTypeNames()[0]).getFeatures();
Encoder encoder = new Encoder(new KMLConfiguration());
encoder.setIndenting(true);
ByteArrayOutputStream os = new ByteArrayOutputStream();
try {
  encoder.encode(features, KML.kml, os);
  String out = os.toString().replaceAll("kml:", "");
  System.out.println(out);
} catch (Exception e) {
  e.printStackTrace();
}

The above code DataStore in = DataStoreFinder.getDataStore(params); is not able to find dataStore and is always returning null.


texts: Detecting Join programmatically using ArcPy?

texts:
I have some Python code that is launched from within an ArcMap project. Any joins that the user may have created in the project must be removed in order for my code to run. Unfortunately, the code that removes a join… arcpy.RemoveJoin_management("layer1", "layer2")… also breaks some of the layer properties that are critical to my application (highlighted fields, read-only fields, etc).

If joins are removed by right-clicking the layer in ArcMap and choosing “Remove Joins” the layer properties are left intact.

If I can detect that a join exists from within my code, I will simply exit the code and display a message that the user must manually remove their joins before attempting to run the code. So… Can a Join be detected programmatically?


texts: Trying to extract xyz from multipolygon shapefile

texts:
I want to extract the x, y and z from shapefile polygon vertices where the geometry is multipolygon.  Once extracted would like to have them in a geopandas data frame. I have seen a couple of approaches however the problem seems to be that when the shapefile is imported I can't get it to a feature that has a Z attribute.  For example, the following works to get the x and y only, using a shapefile that was generated from ArcGIS Pro by converting the multipolygon format to its single parts.
import geopandas as gpd

gdf = gpd.read_file("C:/FG3D_2019/code/Collada_Geospatial/Boston_Mesh_Evaluation/Boston_Bldg1_SingleParts.shp")
print('gdf type', type(gdf))
geom_uu = gdf.unary_union
# Extract the point values that define the perimeter of the polygon
x, y = geom_uu.exterior.coords.xy

What I would like is to get the z coordinates from the geometry field so essentially;
x, y, z = geom_uu.exterior.coords.xyz.
Has anyone had luck with this?


texts:
I have MultiLineString geometry data from an ESRI shapefile.  I also have a polygon (region) that I have created in QGIS using the drawing tools.

When I attempt to find the intersection between the two features using the Geoprocessing tools I end up with an empty set.

If I use the Spatial Query plugin, the features I want are correctly identified but when I try to create a layer I receive the following error: "The query from "YourLayer" using "FID" is not possible".


texts: Is ArcObjects SDK for .NET available with ArcGIS for Home Use or free trial license?

texts:
I currently have a 6 month free trial of ArcGIS Desktop 10.1 that I got online with an authorization number from an Esri Press book and I've read that ArcObjects is bundled with the Desktop license but I don't know if it applies to what I have. I'm not sure that I have the proper installation media with a free trial.

Would the only other way be to gain Customer Service access rights through a Global Account?


texts: How to replace sections of raster with data from another

texts:
I am trying to replace an attribute of one raster dataset, which generalizes the data over that area, with a collection of attributes from another raster which cover the same area, but provide more specific information.

Is there a way to clip out one attribute of a raster, and replace that section with the attributes of another?


texts: non-existing attributes in getfeatureinfo when attribute value is empty

texts:
I am trying to build a table in javascript from getfeatureinfo. I do this by this code:

function showInfo(evt) {
        var f = 0;
        var temstr = "";
        if(evt.features.length != 0){
            for (f = 0; f <= evt.features.length - 1; f++) {
                if(evt.features[f].gml.featureType == "amice_om_herk"){
                    temstr = temstr + "" + "Herkomstlocatie" + "";
                    temstr = temstr + "
";
                    temstr += "
";
                    temstr = temstr += "ID	" +  "Houseletter	" +"City	" +"
";
                    for ( var key in evt.features[f].attributes) {
                        switch (key){
                        case 'afvlstrnum':
                            temstr += ""
                            + evt.features[f].attributes[key] + "	";
                            break;
                        case 'numtv_herk':
                            if (evt.features[f].attributes[key] != null){
                            temstr += ""
                            + evt.features[f].attributes[key] + "	";
                            break;
                            }
                            else {
                                temstr += "	";
                                break;
                            }
                        case 'plts_herk':
                            if (evt.features[f].attributes[key] !== ""){
                            temstr += ""
                            + evt.features[f].attributes[key] + "	";
                            break;
                            }

However, some attributes contain no value. In that case, the getfeatureinfo doesnt recognize the attribute name. So, when i build my table with  and one of the cases for that particular cell is without a value, it skips the case, moves on to the next one and fills the cell with attribute value which should be in the next column. An image of what goes wrong:

Is there a possibility fill  with something like 'nodata' or fill it with nothing? I tried to make an if statement (as seen in the code), but since it never comes into the case when, for that cell, the attribute name doesnt exist, this wont work.


texts: How to expand the raster's extents and add null value into the expansion part?

texts:
I have run a hydraulic model and I got a result of depth raster of the sediment (the brown one). I want to sum the sediment raster with the topography raster (the black one) so I will have the elevation of riverbed after the sediment depth is applied.

My goal is to have a cross-profile of the river after the sediment depth is applied onto the initial riverbed elevation. I used to use Stack Profile tool, and use the help of a line shapefile to sign the location of cross-profile. The problem is that the sediment raster's extent is not large enough to cope with the topography raster so when I created a cross-profile of them, both raster did not align correctly.

What steps should I take if I want to expand the sediment raster's extent and put null values in the expansion cells?

Do you have any idea of what I should do in order to get both cross profiles of sediment raster and topography raster on exact place (for example: the pink line)?


texts:
I have a shapefile containing civic numbers. For each number there are string values of street name and civic number. During the editing it happens that I make typing mistakes and at the end of the work I find some duplicate points, i.e. two or more civic numbers on the same street.

How can I easily find these duplicates in order to export them to another shape using QGIS?


texts: Maximum size on the bounding box with ST_MakeEnvelope and && and geography column in PostGIS?

texts:
I've done quite a few google searches, but haven't managed to find anything that has helped me directly. I've got a few features in a PostGIS table with a geography column. The features are focused around (latitude, longitude) 0, 100 and 1, 101 (those are literally the points for any POINTS or LINESTRINGS in that column).

I'm making a query like this:

SELECT GUID, ATTRS, ST_AsGML(GEOG), ACTIVE, CREATED, FTIME
FROM FEATURES
WHERE GEOG && ST_MakeEnvelope(99, -3, 103, 3, 4326)
 AND ACTIVE = TRUE ORDER BY CREATED DESC LIMIT 100

And it returns all the data

If I change the query to be much larger, then nothing comes back:

SELECT GUID, ATTRS, ST_AsGML(GEOG), ACTIVE, CREATED, FTIME
FROM FEATURES
WHERE GEOG && ST_MakeEnvelope(-75, -80, 110, 3, 4326)
  AND ACTIVE = TRUE ORDER BY CREATED DESC LIMIT 100

So is there a limit I'm hitting somewhere here I don't know about it?

If I change the query to this:

SELECT GUID, ATTRS, ST_AsGML(GEOG), ACTIVE, CREATED, FTIME
FROM FEATURES
WHERE GEOG && ST_MakeEnvelope(-70, -80, 109, 3, 4326)
  AND ACTIVE = TRUE ORDER BY CREATED DESC LIMIT 100

I still get all my data back!

EDIT

Added the ST_AsEWKT(GEOG) as requested:

select ST_AsEWKT(GEOG) from features;
             st_asewkt
-----------------------------------
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;LINESTRING(100 0,101 1)
 SRID=4326;POINT(100 0)
 SRID=4326;POINT(100 0)


texts: TypeError: 'BilFile' object does not support item assignment

texts:
I want to use Python and Spectral library to open a hyperspectral image as ENVI file (.hdr + .dat), rewrite some bands and then save the file in the same format. I have problem with this:

  "TypeError: 'BipFile' object does not support item assignment"

I've got this response when I want to rewrite the band or change the file. Do you know of some solution or advice which could help?

#open file
src = envi.open("C:\soubory_pix\SASI_1.pix.hdr", "C:\soubory_pix\SASI_1.pix")
src_band = src
#rewrite band
src_band[:, 0:319, 96] = src[:, 0:319, 96] * koef_men

#Result
Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: 'BipFile' object does not support item assignment


texts: Creating symbols with multiple colors each from attributes

texts:
I use ArcMap 10.3 to create a query layer of points.
The layer is used in a JavaScript API map as a ArcGISDynamicMapServiceLayer.

Each point has two attributes: point_color and frame_color.
I have 10 different colors for point_color and 10 different colors in frame_color.

The symbol I want is an arrow colored by point_color and a frame colored by frame_color.

I can manually create the symbols as a character marker symbol but then I need to create 10*10 different symbols that deffer only by the colors.

I tried to use LayerDrawingOptions to send a json renderer from the client to the server but there it supports only marker and picture symbols.

Is there a way to create the layer symbology without creating 10*10 different symbols? all symbols are identical except for the point color and frame color.

(I know that I can use a feature layer and render the layer in the client side, but the layer is VERY large and I have performance issues so I want to render the layer is the server side)

Here is a sample of a symbol with red as frame_color and green as point_color:


texts: Replicas: Synchronize changes for individual feature class (not entire replica)

texts:
I have two ArcSDE databases, running with a two-way replica.One serves our interactive mapping website which gives certain users the ability to update, add points, edit attributes, etc. Then we have our master which is updated frequently by power users.

There are only a few feature classes that are currently being edited, and  I would like to synchronize these updates ONLY if these feature classes have been edited.

I've developed a python script using the compare features tool, which compares the feature datasets with the master. I'd like to only synchronize changes from these three layers, update the master (so child to parent) and then have the master sync any new changes to the edit database (parent to child).

Is there any way to synchronize changes and schema for individual datasets or is there a better way to go about this?

import arcpy
aburg_master_DBO_TownOwnedJan2013 = "Database Connections\\Connection to       gisserver.sde\\aburg_master.DBO.Geocortex\\aburg_master.DBO.TownOwnedJan2013"
Aburg_Edit_DBO_TownOwnedJan2013 = "Database Connections\\Connection to gisserver(2).sde\\Aburg_Edit.DBO.Geocortex\\Aburg_Edit.DBO.TownOwnedJan2013"
compare_txt = "C:\\Test\\compare6.txt"

# Process: Feature Compare
result = arcpy.FeatureCompare_management(aburg_master_DBO_TownOwnedJan2013,   Aburg_Edit_DBO_TownOwnedJan2013, "GlobalID", "ATTRIBUTES_ONLY", "IGNORE_M;IGNORE_Z", "0.001 Unknown", "0", "0", "", "OBJECTID;created_user;created_date;last_edited_user;last_edited_date", "NO_CONTINUE_COMPARE", compare_txt)
changes = result.getOutput(1)
print "Comparision Result is: " + changes
change = str(changes)
if change == 'false':
    print "LOL"
else:
    print "Nope"


texts:
In QGIS I'm trying to label the neighbouring features of a certain polygon but not getting much success. All the features are inside the same layer named MG_Municipios_2020. I'm using the Rule-Based Labeling with this expression to verify which geometries intersects with the one I want so that I can lable them:
intersects (
      $geometry, geometry (
             get_feature(
                 'MG_Municipios_2020', 'MN_MUN','Coronel Fabriciano'
             )
      )
)

Can somebody help me to understand why is it not working?


texts: Outputting geoJSON with arcpy.FeaturesToJSON_conversion tool of ArcGIS Pro?

texts:
I'm using ArcGIS Pro 2.0 and I've updated my system path to point to the version of Python (3.5.3) which shipped with ArcGIS Pro.

FeaturesToJSON_conversion is well documented as including a parameter for exporting to geoJSON in the new version.  The desktop version of ArcPy does not include this parameter which is why I upgraded.

I receive a 999999 error when specifying geoJSON='GEOJSON'. Just to be sure I'm using the right version of ArcPy I confirmed that it can take the parameter geoJSON='NO_GEOJSON'

I know how to do this with ogr2ogr -f GeoJSON but I was hoping to stay in ArcPy this week.

This gives error 999999:

arcpy.FeaturesToJSON_conversion("lyr", r"G:\GIS_Projects\Backups\json\myjsonfeatures.json", geoJSON='GEOJSON')

but this works fine:

arcpy.FeaturesToJSON_conversion("lyr", r"G:\GIS_Projects\Backups\json\myjsonfeatures.json", geoJSON='NO_GEOJSON')


texts: Create a combobox in qgis plugin to input text (string)

texts:
I'm developing a plugin in QGIS, with a dialog box where I want to insert text (string) and send to web service soap.

I don't found much on the internet.

Except this code that i have adapted:

qid = QInputDialog()
canvas = iface.mapCanvas()
input, ok = QInputDialog.getText( qid, "Enter Date", "Enter text", QLineEdit.Normal, "Date" + "," + "Text")
if ok:
    Date = input.split( "," )[ 0 ]
    print Date
    Name = input.split( "," )[ 1 ]
    print Name

But this code don't work in the combobox that i want.

I create this ComboBox in QCreator:

Can someone help me to adapt the code to this combobox?


texts: How to remove part of a string value in the attribute table of a shapefile in PyQGIS

texts: I have a number of shapefiles with a field in the attribute table containing many different string values similar to this one: 'Waypoints_17-JUL-18Timothy waypoints.shp Waypoints_17-JUL-18Timothy waypoints.' The only difference in the values is the name. I would like to find a way in PyQGIS to remove everything after the first name, so before 'waypoints.shp...' and iterate through each of the fields to repeat this. Is it possible to do this in PyQGIS? I am new to PyQGIS so am unsure where to start but it should be a lot quicker than doing manually. I am using QGIS 3.0.0


texts: Execute a list of SQL statements (ArcPy)

texts:
I'm just starting out with ArcPy and ArcSDESQLExecute.

I'm trying to use the 'Execute a list of SQL statements' ArcPy script from http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//002z00000021000000.htm

import arcpy
from arcpy import env
import sys

try:
# Make data path relative
#
env.workspace = sys.path[0]

# Two ways to create the object, which also creates the connection to ArcSDE.
#   Using the first method, pass a set of strings containing the connection properties:
#   <serverName>,<portNumber>,<version>,<userName>,<password>
#   sdeConn = arcpy.ArcSDESQLExecute("gpserver3","5151","#","toolbox","toolbox")
#   Using the second method pass the path to a valid ArcSDE connection file
#
sdeConn = arcpy.ArcSDESQLExecute(r'Database Connections\GIS_Connection.sde')

# Get the SQL statements, separated by ; from a text string.
#
SQLStatement = arcpy.GetParameterAsText(0)
SQLStatementList = SQLStatement.split(";")

print "+++++++++++++++++++++++++++++++++++++++++++++\n"

# For each SQL statement passed in, execute it.
#
for sql in SQLStatementList:
    print "Execute SQL Statement: " + sql
    try:
        # Pass the SQL statement to the database.
        #
        sdeReturn = sdeConn.execute(sql)
    except Exception, ErrorDesc:
        print ErrorDesc
        sdeReturn = False

    # If the return value is a list (a list of lists), display each list as a row from the
    #   table being queried.
    if isinstance(sdeReturn, list):
        print "Number of rows returned by query: " + len(sdeReturn), "rows"
        for row in sdeReturn:
            print row
        print "+++++++++++++++++++++++++++++++++++++++++++++\n"
    else:
        # If the return value was not a list, the statement was most likely a DDL statment.
        #   Check its status.
        if sdeReturn == True:
            print "SQL statement: " + sql + " ran sucessfully."
            print "+++++++++++++++++++++++++++++++++++++++++++++\n"
        else:
            print "SQL statement: " + sql + " FAILED."
            print "+++++++++++++++++++++++++++++++++++++++++++++\n"

except Exception, ErrorDesc:
    print Exception, ErrorDesc
except:
    print "Problem executing SQL."

I'm getting this error:
" 'ascii' codec can't encode character u'\udf9c' in position 52: ordinal not in range(128)"

What have I done wrong? I don't know if it's a problem with the ascii formatting in the script, or if the SQL statement needs to be put in the script somewhere, or if it's something else?

ArcGIS/SDE 10.3.1


texts: Labels are diappearing while zooming in or out

texts:
We are using Geoserver 2.10.0. While displaying the labels in the maps, the labels which are there on the border of the tiles(4*4) are not getting displayed properly. the labels can be seen after further zooming in or out(when they fall in the center of the tile and not in the border). in this image you can see the labels are being cut

On zooming further, the labels come properly, as expected

In the first image(1.png), the labels are cut. But in the second image(2.png), we see on zooming further, the labels are coming as expected.

Please help.

Can someone please suggest a solution without using GeoWebCache.


texts: Leafletjs popups in iframe not working in IE11

texts:
This is my leafletjs iframe in a wordpress blog: http://www.kemeling-consulting.com/cresh-map-viewer/

The iframe is working fine in firefox and chrome, in IE11 the popups are not working. They are getfeatureinfo requests from the geoserver WMS.

The website also works fine (popups work ok) in IE11 directly: http://185.19.29.113:8080/cresh_map/mapla.html

I added this in the header.php file

Does anybody know what else I can try?


texts: How to use ArcPy to list feature class and number of rows per feature class?

texts:
I would like to iterate a geodatabse to get the name of each feature class/data set in it, in addition to the number of rows in the feature class. I would also like to identify whether the dataset that is being printed is a feature class or feature dataset.

The code that I have now prints each feature class and feature data set; but does not denote whether if the object printed is a feature class or feature data set.  How would I print the number of rows per feature class in addition to denote whether if it is a feature class or feature dataset?

I am using ArcGIS Desktop 10.0 and a direct connect via spatial database.

import arcpy
arcpy.env.workspace='myworkspace'

    datasetList = arcpy.ListDatasets("*", "Feature")
    for dataset in datasetList:
        print dataset
        fcList = arcpy.ListFeatureClasses("*","",dataset)
        fcList.sort()
        for fc in fcList:
            print fc
        for fc in fcList:
            arcpy.GetCount_management(fc)

Sample Output:

SDEDEV.SDE.Boundaries
SDEDEV.SDE.CTYA_GLEN
SDEDEV.SDE.ZipCode
SDEDEV.SDE.Streets
SDEDEV.SDE.CAMS_Streets_2012
SDEDEV.SDE.ROW
SDEDEV.SDE.Streets_COG


texts: Using raster2pgsql to display data in QGIS

texts:
I am quite new to this process (using raster data in PostGIS) and tried to find some tips/solutions on other threads but I didn't manage to get things working.

I ingested some Landsat data into my database using:

raster2pgsql -I -C -s 28356 LS7_20150425_B50.tif public.LS7_20150425_B50 | psql –d postgres –U postgres

And everything is working fine : I see my new table LS7_20150425_B50 with two columns (rid and rast) and I am able to compute some stats on my raster using ST_SummaryStats().

However, I'd like to display my raster in QGIS (2.12.0 Lyon) using the DB Manager. However, when I am adding the table to the canvas I am getting an error message:

PG: dbname=postgres host=localhost user=postgres port=9999 mode=2 schema=public column=rast table=ls7_etm_nbar_p54_ganbar01-002_089_078_20150425_b40 is not a valid layer.

Do you have any idea?


texts: How to transform roadline uid to point geometry?

texts:
I have a road line geom and roadline uid attribute. How to transform selected roadline uid to selected points geoms?

ArcGIS has spatial adjustment tool, does QGIS have a similar tool?

A spatial join doesn't solve the problem.


texts:
I've set up a tile rendering system based on Mapnik + PostGIS. When I open the data with QGiIS directly from PostGIS it looks fine; but when I attempt to generate tiles with Mapnik they look stretched (notice that this issue only appears when the rendered tiles are on another horizontal line).

Tiles viewed with OpenLayers:

Tiles viewed with QGIS:

Raw data from PostGIS viewed with QGIS:

Can you please guide me trough what things should I check in order to fix this issue?


texts: Extract all pixels from a finer DEM into a coarser fishnet in each grid box in arcgis

texts:
I have a raster DEM of 90m resolution (srtm90). I have a fishnet(polygon) for the same region of size 25*25 sq. km. For each grid in fish net, I want to store all the pixel values of raster within the same lat-lon range of fishnet grid box. The ultimate aim is to create elevation bands of different range say 0-1000m, 1000-1500m, 1500-2000m and so on with their respective areas inside each fishnet grid box.

The following picture states the problem with a table asking for each elevation pixel with its respective area. Though area can be same for each pixel in a DEM, all pixel values are required in each grid box. These area are later used to determine the area fraction of each elevation bands inside a grid.

Raster to polygon followed by intersection with fishnet was my initial attempt but the file was so big that raster to polygon conversion was never going to stop. I dont think this method is appropriate. Hence I am bound to question here.

Arcgis is preferred.


texts: QGIS - Vector layer rendering - dissolving objects

texts:
I have a vector line layer for a stream network.  I'm styling the layer with a 100 map unit line width with 50% opacity.  Opacity is set to 100% within the symbol class and opacity is set in the Layer Rendering section.  The lines are not continuous and there are overlaps along streamlines and where tributaries join.
At these overlaps the opaque lines are darkened. The effect is identical if using geometry generator ro buffer by 50m (rather than a simple 100m line width).  A workaround would be create a seperate buffer layer and dissolve the features, but I want to avoid creating unnecessary additional layers that would lose their link to the underlying stream network.
I've cycled through the options uder Feature blending mode but cant find an option that hides darkening over overlaps.
Is there a setting in layer rendering that will allow overlaps to be ignored within that layer?


texts:
I have 2 huge files of points (hundreds of thousands of points).  I need to determine a way to loop through all the points in one file and find all the points in the second file that are only a certain distance away.  So like within a distance outside the buffer.  Or on the outer edge of the buffer.

I am thinking of 2 buffers, one larger buffer and then one smaller buffer and then somehow subtracting the values.  No idea how this could be accomplished though.

The files are huge so I am also concerned in doing this the quickest way possible.

Any ideas?


texts: Linking PDF file to a polygon shape with relative paths?

texts:
I´m looking for a solution to connect a specific PDF file to each polygon in my shape and to create an action to open the PDF. Relative paths are required.

I am aware many people already asked this kind of question, but they´ve been
asked and answered years ago.

I´m using QGIS Madeira 3.4.11 (long term release) and Windows 10 (64 bit).

My file structure is like this:

Gemeinde.qgs

-Daten

--Oekokonto

---Datenblaetter

----10.pdf

With absolute paths I did it like this:

action type: open

action: [%Datenblatt%]

The absolute filepath is in the column "Datenblatt" in my attribute table.

Example: X:/Projekte/QGIS/Daten/Oekokonto/Datenblaetter/10.pdf

What do I have to change to have relative instead of absolute paths?

Edit:

Attribute table:


texts: Why is my arcpy.Polygon changing how many vertices it has?

texts:
Using a Python script, we are pulling a list of polygon vertices from a web service to build a polygon and insert it into our SDE.

The problem steps are happening when we're making an arcpy.Polygon geometry object out of an arcpy.Array of arcpy.Points. There are 95 points listed in the JSON, 95 points in the polygon_array... and then not 95 points in the resulting arcpy.Polygon.

If it's created with no spatial reference, the polygon pointCount has 11 points
If it's created with WGS84 (which matches the lat/long projection of the source data), the polygon pointCount has 1126 (!) points

Why does the resulting polygon have such a different number of vertices (either less or more) than it's being handed? Is it possible for misordered vertices that would cause geometry problems to be "weeded out" by the arcpy.Polygon instantiation — but if so, I'm even more baffled by the WGS84 polygon having over a thousand points!

Code below:

# add points to array
polygon_array = arcpy.Array()
for point in polygon_data:
    polygon_array.add(
        arcpy.Point(
            polygon_data[point]['Value']['Longitude'],
            polygon_data[point]['Value']['Latitude']
        )
    )

# print the length of the array of points
print(polygon_array.count) # prints 95

sr = arcpy.SpatialReference(4326)
polygonObj = arcpy.Polygon(polygon_array, sr)

# Print the amount of points in our polygon
print(polygonObj.pointCount) # prints 1126

(This is running with ArcGIS Desktop 10.2.1.)


texts: Filling gaps with two conditions, closest existing data and same name

texts:
I'm working with a traffic dataset, with four main fields: name of the streets, intensity and lat-lon coordinates. So it is something like this (totally made up data)

street   IMD     lat     lon

a         27     41.423  2.121

a         NA     41.503  2.088

a         35     41.214  2.102

b         NA     40.231  2.484

b         13     40.235  2.471

b         29     40.254  2.814

I need to fill the NA data with the closest imd data, and also it must be the same street, so two conditions. And if there is a data already, don't do anything. I'm trying with "for" and "if" lines, but I'm totally lost. Of course we are talking about thousands of lines.

I've been looking for an answer for this but I've not found it... I'm still a beginner in R, but I'm doing my best.


texts:
I am using QGIS and the SCP plug-in to classify a sentinel 2 image.

As for the reference shapefile my point layer, is a vector point layer containing the names of crops, which were used in creating ROIs.

My question is:
Do I need to alter this crop point file file further in order to generate a valid accuracy assessment? I am aware it needs to be converted into either a shapefile or raster, but
Do I need to attach the pixel values from the ground truth image to the crops, in order to use this as a reference layer?

I have computed accuracy assessments already but they are between 0.02 and 4%.


texts: How to simulate events from Openlayers 2?

texts:
I am working on a GIS API that uses OpenLayers. I am testing it using QUnit.

The problem is I need to simulate clicks and trigger some events to be able to test the whole code.

For example, I have a feature to draw a polygon and I want to just trigger the "featureadded" event. Is it possible to trigger OpenLayers events without actually working graphically?

Thanks


texts: Error creating a spatial database. ERROR

texts:
I have been trying to get Postgis set up for a while now, I am stuck on one error and have not found a solutions yet. This is on Ubuntu 14.04.1.

I am trying to run:

CREATE EXTENSION postgis

But am getting this error:

ERROR: could not load library "/usr/lib/postgresql/9.1/lib/rtpostgis-2.0.so": libgdal1.7.0.so

I am taking it's saying I need to get libgdal installed but every thing I have found so far doesn't seem to work.


texts: Taking 2 overlapping spatial polygon layers and creating 1 with features from one replacing those of other layer?

texts:
I have two sf multipolygon layers - one is a 5km buffer along the Mexico-US border and the other is an administrative 2-level layer of Mexico. I am trying to figure out how to layer these on top of one another and have the final product be one layer that can be visualized like this:

I have tried spatially joining with st_join, however this results in either just the border polygon with attributes from the admin 2 layer (listing buffer layer first) or in the admin 2 layer with buffer layer attributes. I am looking to get one layer with the polygons from both and attributes in overlapping areas to be defined by the buffer layer.


texts: ArcGIS JavaScript API layerDrawingOptions seem restricted to 1 option, and it doesn't have any effect

texts:
I have two problems with ArcGIS Server's (10.0) dynamic layers when calling Layer.setLayerDrawingOptions from v3.2 of the JavaScript API.

1) It doesn't work. I set the layerDrawingOptions array for a dynamic layer and the resulting image is rendered the same as if I hadn't changed anything. I have verified that the layer ID is correct, and also tried several approaches to defining the LayerDrawingOptions object.

When I inspect the request that goes to my map server it contains the following dynamicLayers property:

[{"id":0,"source":{"type":"mapLayer","mapLayerId":0},"definitionExpression":"ASSET_ID IS NULL"},{"id":1,"source":{"type":"mapLayer","mapLayerId":1},"definitionExpression":"ASSET_ID IN (27)","drawingInfo":{"renderer":{"type":"simple","symbol":{"color":[255,0,0,255],"outline":{"color":[0,0,0,255],"width":1,"type":"esriSLS","style":"esriSLSSolid"},"type":"esriSFS","style":"esriSFSSolid"}},"transparency":1,"scaleSymbols":false,"showLabels":false}},{"id":2,"source":{"type":"mapLayer","mapLayerId":2},"definitionExpression":"ASSET_ID IS NULL"},{"id":3,"source":{"type":"mapLayer","mapLayerId":3},"definitionExpression":"ASSET_ID IS NULL"}]

The feature with ASSET_ID = 27 is rendered but in its default style, however "color":[255,0,0,255] says it should be solid red. The map service referenced by this layer only contains a single polygon layer and that is ID 1. I can't see any way that I am accidentally attempting to style the wrong layer.

2) If I define multiple layerDrawingOptions in the array, one for each layer within the map service, I am forced to provide a proxy URL (even though the resulting request is nowhere near the 2000 character limit). When what should be the map request goes to the map server via the proxy, with multiple layerDrawingOptions, the map server returns a JSON string that describes an image (height, width, bbox etc.) but has an empty href ("").

Although there are no errors in my AGS log I assume this means I'm not allowed to set multiple layerDrawingOptions in a single layer. Is this actually true?

So if anyone can provide thoughts or experience on 1) why layerDrawingOptions seem to have no effect, or 2) why I am seemingly not allowed to set more than one option for a map service, I'd really love to hear it/them.


texts:
I have this image that I want to convert it to GeoJSON.

And I have this World File:

0.0388048118
0.0
0.0
-0.0387897595
697054.11331
4458020.56656

When I try to run the gdal_polygonize command from terminal, it pixelizes the json, but If I try to run it without the World File it works very well:

gdal_polygonize.py -f "GeoJSON" image.png out.json

What can I do to make it work with the World File?
I need the coordinates, so I can't run it without that file...
Is there any alternative to gdal_polygonize?

EDIT

This is what I should get. This GeoJSON has been created without the World File:

This is what I get with the World File:


texts: How to toggle layer visibility in QGIS using Python?

texts: How can I toggle layer visibility in QGIS via Python API?


texts: special python library needed for spatialite?

texts: It appears that pyspatialite is very outdated.  Do I need a special library like pyspatialite in order to use the spatial functions of spatialite?  Or can I just connect to the database via the normal python driver and execute the spatial functions?  If so, is there a more updated library that I am not aware of?


texts: does the census bureau provide clipped / water-bounded files as nice as certain cities?

texts:
this question says they do but when i compare what i can find in the cartographic boundary files section of the census bureau's gis page, they don't look as good as certain shapefiles provided by the city government.  am i missing something?  or do i have to scrounge through each city's website to get this level of detail (assuming it's there at all)

comparison photos of new york city and reproducible R code below.

library(maptools)

tf <- tempfile() ; td <- tempdir()

# from http://www.census.gov/geo/maps-data/data/cbf/cbf_cousub.html

download.file( "http://www2.census.gov/geo/tiger/GENZ2013/cb_2013_36_cousub_500k.zip" , tf )
z <- unzip( tf , exdir = td )
ny.shp <- readShapePoly( z[ grep( 'shp$' , z ) ] )
nyc.shp <- subset( ny.shp , as.numeric( as.character( COUNTYFP ) ) %in% c( 5 , 47 , 61 , 81 , 85 ) )
# shows some water.
plot( nyc.shp )

# from http://www.nyc.gov/html/dcp/html/bytes/districts_download_metadata.shtml

download.file( "http://www.nyc.gov/html/dcp/download/bytes/nybb_14c.zip" , tf )
z <- unzip( tf , exdir = td )
nyc.shp <- readShapePoly( z[ grep( 'shp$' , z ) ] )
# very awesome.  does not show water.
plot( nyc.shp )


texts: ogr2ogr{gdalUtils} correct syntax to set Ground Control Points for vector georeferencing

texts:
I would like to know correct syntax of R {gdalUtils} package (version 2.0.1.7) to georeference Shapefiles using corner points as GCPs. (These Shapefiles have y=height and x=width information, without datum / projection.)

On OSGeo4W shell, ogr2ogr works fine and it returns georeferenced Shapefile.

  ogr2ogr -gcp 0 0 40 40 -gcp 0.263693 0 145 40 -gcp 0.263693 -0.237619 145 -50 -gcp 0 -0.237619 40 -50 -f "ESRI shapefile" Myreferenced.shp My.shp

Trying to do the same by gdalUtils::ogr2ogr(),

library(gralUtils)
src_datasource_name # input Shapefile
dst_datasource_name # output Shapefile (georeferenced)

# ogrinfo() works fine
ogrinfo(datasource_name = src_datasource_name, al = TRUE, so = TRUE)

# ogr2ogr() fails
ogr2ogr(src_datasource_name= src_datasource_name,
        dst_datasource_name= dst_datasource_name,
        f= "ESRI Shapefile",
        gcp = list(1:4,
                   c(0, 0, 40, 40),
                   c(0.263693, 0, 145, 40),
                   c(0.263693, -0.237619, 145, -50),
                   c( 0, -0.237619, 40, -50)
                  )
       )

What I am struggling in the ogr2ogr() syntax is how to provide multiple GCP points.
As you see above, I tried to give GCPs as a list, which did not work.

Does anybody know how to set GCPs with 4 corner points in gdalUtils? Unfortunately use cases of gdalUtils are relatively scarce on internet.

About my Shapefile Data

They were originally hand-drawn by MS-Powerpoint and saved as .eps file. It has bounding box with known LatLong corner points (UL, UR, LR, and LL). The .eps file was then converted to .shp file.
Conversion from .eps to .shp changes everything into Line feature (Multilines Shapefile).
This Shapefile does not have projection information but its dimension can be read by ogrinfo. It is (0, 0) at the UL corner, with height (eg y= -0. 238) and width (x= 0.263). Measured unit is meter. These y and x values can vary among drawings.

Error message from above ogr2ogr() is

[1] "ERROR 1: Failed to compute GCP transform: Not enough points available"

However increasing points does not seem a solution, as ogr2ogr on OSGeo4W shell works fine on the same data and GCP sets. This error is due to my wrong code, but I do not know how to fix it.

[EDIT] by QGIS Raster georeferencer, .points would look like this image.

mapX,mapY,pixelX,pixelY,enable
40,40,0,0,1
145,40,0.26300000000000001,0,1
40,-50,0,-0.23799999999999999,1
145,-50,0.26300000000000001,-0.23799999999999999,1

[EDIT2] Another way of applying GCPs:

ogr2ogr(src_datasource_name= src_datasource_name,
        dst_datasource_name= dst_datasource_name,
        f= "ESRI Shapefile",
        gcp = c(0,                0,  40,  40),
        gcp = c(0.263693,         0, 140,  40),
        gcp = c(0.263693, -0.237619, 140, -50),
        gcp = c(0,        -0.237619,  40, -50)
        )

This way R returned following error

  Error in ogr2ogr(src_datasource_name = src_datasource_name, dst_datasource_name = dst_datasource_name,  :
    formal argument "gcp" matched by multiple actual arguments

I don't know what it means...

[EDIT3]
Suggested approach to put variables into -optfile.

ogr2ogr(src_datasource_name= src_datasource_name,
        dst_datasource_name= dst_datasource_name,
        f= "ESRI Shapefile",
        optfile = "C:/temp/my_gcp.txt"
        )

Below is my tested formats (my_gcp.txt) yet successful...

# trial 1
"-gcp 0 0 40 40 -gcp 0.263693 0 145 40 -gcp 0.263693 -0.237619 145 -50 -gcp 0 -0.237619 40 -50"
# trial 2
"gcp = c(0, 0, 40, 40,), gcp = c(0.263693, 0, 145, 40), gcp = c(0.263693, -0.237619, 145, -50), gcp = c(0, -0.237619, 40, -50)"

Error message was:

Error in ogr2ogr(src_datasource_name = src_datasource_name, dst_datasource_name = dst_datasource_name,  :
unused argument (optfile = "C:/temp/my_gcp.txt")


texts: Using URL parameters in ArcGIS Server Viewer for Flex?

texts:
Can you use URL parameters in ArcGIS Server 10 Viewer for Flex to turn on layers?

I have had success using URL parameters to go to specific coordinates.


texts: Shapefile to PostGIS via shp2pgsql.exe failed

texts:
I am trying to import a shapefile to PostGIS via shp2pgsql.exe.
I am not successful. I am only get the following message:

Importing with configuration: tune_aoi, public, geom, C:\Temp\TUNE_AOI.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=0
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Shapefile import failed.

Here's a link to the dataset: https://www.dropbox.com/s/325qlkeg3pez100/Temp.zip?dl=0

What am I doing wrong?


texts: Import non spatial data from MySQL DB into QGIS using queries

texts:
I want to import data from MySQL DB.
MySQL schema has a lot of tables with a huge amount of data. This is non spatial data. I want to import objects as vector layers using sql query which is selecting and joining only what I need.
Is it possible to do?

I`ve tried to use Layer->Add Vector Layer ...> MySQL database. But QGIS get crashed in this case (Import MySQL data to QGIS with a plugin).
Indeed, my db has not vector data, it has only tables with longitude, latitude columns.

I`ve tried to use DB Manager plugin but its help says that only four predefined dbms backends are supported (db_manager github) and no MySQL DBMS listed there.


texts:
Usually when I create maps with a set of numbered points (e.g. numbered locations on the map which link across to a separate table describing these points) I try to number them in a rough "reading order". I.e. in more-or-less horizontal strips from left to right, beginning at the top of the page.

I've started wondering lately if this is the best approach, but am unable to find any guidelines or research on what is the most ideal numbering sequence for maps.

So now I'm seeking any best-practice guides or publications discussing this. (Obviously the main guiding factor here is aiding map readers to quickly locate a specific numbered point on the map itself!)


texts: Arguments not working in cartoDB createLayer method

texts:
I added a cartoDB layer to an existing map using createLayer method and it loaded properly. Thereafter, created a variable params and parsed the arguments as sees on the code snippet below. Unfortunately I am not able to see the title, shareable button and the search bar on my map. What could be wrong?

var params = {
  search: true,
  cartodb_logo: false,
  mobile_layout: true,
  shareable: true,
  title: true
}

cartodb.createLayer(map, jsonUrl, params).addTo(map)
.on('done', function(layer){
  layer.setZIndex(1000); baseLayer
})


texts: Is there a one step operation to select all the records in a QGIS layer?

texts:
Is there a one step operation to select all the records in a layer or dataset in QGIS? I can't find any reference to a single button, menu command or keyboard shortcut to achieve this.

I am aware of workarounds such as the following:

Selecting one record in the attribute table, inverting the selection and then reselecting the now unselected record
Running an advanced select with a static 'True' condition such as 1=1.
Using a marquee or other polygon select tool over the entire layer extent in the map window

I'm running QGIS 1.7.4 on Windows 7


texts:
GeoServer's web service base URL is http://localhost:8080/geoserver/www/.
How can I change it to be just http://localhost:8080/ in order to link a domain name to my machine?


texts:
We have installed our own OSM overpass server with planet data. It was first working ok, but after enabling data updates it slowed down significantly. CPU and disk IO usage is low. I cannot say for sure that data updates is the cause. Is there something describing how to troubleshoot it? https://wiki.openstreetmap.org/wiki/Overpass_API/Installation
has some info, but not about troubleshooting performance.

Hardware info: m4.2xlarge instance type in AWS with GP2 SSDs. GP2 = general purpose SSDs. 8 vCPUs, 31GB RAM.

M4 instances are the latest generation of General Purpose Instances. This family provides a balance of compute, memory, and network resources, and it is a good choice for many applications.
Features:
2.3 GHz Intel Xeon® E5-2686 v4 (Broadwell) processors or 2.4 GHz Intel Xeon® E5-2676 v3 (Haswell) processors
EBS-optimized by default at no additional cost
Support for Enhanced Networking
Balance of compute, memory, and network resources
Model    vCPU    Mem (GiB)    SSD Storage (GB)    Dedicated EBS Bandwidth (Mbps)
m4.large    2     8           EBS-only            450
m4.xlarge   4    16           EBS-only            750
m4.2xlarge  8    32           EBS-only           1,000


texts: How to Interpolate through ModelBuilder?

texts:
I would like to create an iterated interpolation tool through ModelBuilder. I have an excel file containing different columns and each column represents a data that will be interpolated using Spline Tool. How can I do this through ModelBuilder? I am using ArcGIS 10.

Image contains point data with lat/long and sample# are the columns to be interpolated.


texts: Difference between contrast enhancement options QGIS 3 rasterlayer symbology multiband color

texts:
When symbolizing a multiband color raster layer, in QGIS there is a field called contrast enhancement. In this field I have four options to choose from:
'no enhancement'  'Stretch to MinMax' 'Stretch an Clip to MinMax' 'Clip to MinMax'. What is contrast enhancement good for and what is the difference between the options?


texts: Creating interface for Editing ArcSDE Layers in web service online?

texts:
The initial idea is to make a multiuser web service for editing versions of an ArcSDE geodatabase online (each version for each user. Is it possible to acces through a password and user name?).
Can I create my own personal interface with specífic tools of ArcGIS?

I need to edit 3D vertices of lines, polygons, points, create new features, copy, paste, clip, delete, make selections,... and other simple tools, also edit the attribute tables like in ArcGIS Desktop (but not geoprocessing tools), and I don't know if this is possible.

Can anyone give me an idea of what can I create and the steps I need to study and prepare this idea?


texts: Storing image with rasterio as gtiff modifies image histogram

texts:
I'm trying to store and image as as gtiff, but the resulting image is changed, and I can't explain why.

I'm reading sentinel jp2 band images using rasterio, I rescale these images from [0-10.000] to [0-255].
Four bands are combined to create a [x,y,4] numpy array.

After that, the image is processed keeping the range [0-255] using numpy/opencv.

I want to store the resulting image in a gtiff file, but the image is modified during the process, and the original and resulting histograms are much different.

This is the profile from one of the original sentinel images:

  {'driver': 'JP2OpenJPEG', 'dtype': 'uint16', 'nodata': None, 'width': 10980, 'height': 10980, 'count': 1, 'crs': CRS.from_dict(init='epsg:32630'), 'transform': Affine(10.0, 0.0, 300000.0,
         0.0, -10.0, 4600020.0), 'blockxsize': 1024, 'blockysize': 1024, 'tiled': True}

To store the resulting image, I use this profile:

  {'driver': 'GTiff', 'dtype': 'uint8', 'nodata': 0, 'width': 10980, 'height': 10980, 'count': 4, 'crs': CRS.from_dict(init='epsg:32630'), 'transform': Affine(10.0, 0.0, 300000.0,
         0.0, -10.0, 4600020.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'interleave': 'pixel', 'compress': 'jpeg'}

The code

# Here i'm reading one of the original sentinel images to get it's profile,
# change it, and use it to geolocate the processed image
dataset = rasterio.open('one_of_the_original_sentinel_images')
profile = dataset.profile
dataset.close()

# change original profile to set geotiff driver
profile.update(
    **DefaultGTiffProfile(dtype=rasterio.uint8, count=img.shape[-1],
                          compress='jpeg', interleave='pixel',
                          # photometric='ycbcr', can't use it with 4-band image
                          ))
# store the processed image (img) into a 4-band raster
with rasterio.open(img_filename, 'w', **profile) as dst:
    for ch in range(img.shape[-1]):
        # iterate over channels and write bands
        img_channel = img[:, :, ch]
        dst.write(img_channel, ch + 1)  # rasterio bands are 1-indexed

This is the image histogram before storing as gtiff with rasterio

Image histogram after storing it as gtif with rasterio (optained reading the gtif file directly with opencv.

UPDATE

Sorry, talking about the sentinel jp2, can be confusing mislead away from the main problem. I included it because I think the differences in profiles is the key of the problem.

My problem: I have an image stored in a numpy array, each channel in the range [0-255], and I want to geolocate it and store it in a gtiff file.
To obtain he csr and affine, I use the original sentinel profile because I want the resulting image to have exactly the same geoposition.

The histograms I show above belong the image before and after storing it using rasterio. Besides, if I store the image directly using opencv, the image is right.
So I think the pixel length or data range is not the cause of the problem.

So, which is the right way to get he width, height, affine and csr from a raster and use it to create a new raster with an image of the same size?

new_profile = DefaultGTiffProfile(dtype=rasterio.uint8, count=img.shape[-1], compress='jpeg',
                                      interleave='pixel')

# change the w/h, crs and affine from the sentinel original images
new_profile.update(crs=profile["crs"], transform=profile["transform"],
                   width=profile["width"],
                   height=profile["height"])

with rasterio.open(img_filename, 'w', **new_profile) as dst:
    for ch in range(img.shape[-1]):
        # iterate over channels and write bands
        img_channel = img[:, :, ch]
        dst.write(img_channel, ch + 1)  # rasterio bands are 1-indexed


texts: Selecting features that intersect features of another polygon with spatial index

texts:
Based on the answer of @mgri here: https://gis.stackexchange.com/a/224954/164953 I tried to select the features of a multilines layer intersecting less or more than 2 features of a multipoints layer.
At first sight, I thought that the script did what I wanted but unfortunately, I realized that there were a significant number of features that were (not)/selected whereas they (should)/shouldn't
from qgis.core import *
from qgis.gui import *
from qgis.utils import *

layer1 = QgsProject.instance().mapLayersByName('lines')[0]
layer2 = QgsProject.instance().mapLayersByName('nodes')[0]

index = QgsSpatialIndex() # Spatial index
for ft in layer2.getFeatures():
    index.insertFeature(ft)

selection = [] # This list stores the features which contains at least one point
for feat in layer1.getFeatures():
    inGeom = feat.geometry().buffer(0.1,0.1)
    idsList = index.intersects(inGeom.boundingBox())
    if len(idsList) != 2:
        selection.append(feat)

layer1.selectByIds([k.id() for k in selection])

A small sample of the dataset: https://we.tl/t-caDtSeAe1z


texts: Finding Points where two polylines cross in ArcGIS Runtime .NET SDK?

texts:
I am using ArcGIS SDK for .NET and have a graphic overlay defined on my mapview.  I added two circles to the map (made with polylines) which overlap.  I now need to determine the points at which the two polylines cross.

The GeometryEngine.Intersection() method returns an empty geometry.  The GeometryEngine.Crosses() method returns true.

How can I get the points at which the two polylines cross?

Edit:  Here is some code that I was trying:

' Get reference to overlay for add operation
Dim overlay As GraphicsOverlay = Me.ehMapControl.MapView.GraphicsOverlays.Item("NewLocation")
' Create a graphic object
Dim myGraphic as Graphic = New Graphic
' Create a polyline that is a 15 foot circle around "point"
Dim myEllipseParams As New GeodesicEllipseParameters(point, 15,  LinearUnits.Feet)
myEllipseParams.OutGeometryType = GeometryType.Polyline
myGraphic.Geometry = GeometryEngine.GeodesicEllipse(myEllipseParams)
' Add polyline to graphics overlay
overlay.Graphics.Add(myGraphic)
' Create a second polyline
Dim myGraphic2 as graphic = new Graphic
...
' Perform an intersection
Dim intsersect = GeometryEngine.Intersection(myGraphic, myGraphic2)
' Doesn't work (returns an empty geometry)
Dim isCrossing = GeometryEngine.Crosses(myGraphic, myGraphic2)
' Returns true, so the engine knows the two polylines cross


texts: How should variable data related to a point be stored in a database table?

texts:
I've been having a bit of trouble understanding how a dataset is built/managed for GIS applications.
Over the past few weeks (part time) I've been trying to understand the basics of GIS, and how it works with the likes of GeoServer, QGIS, etc... in terms of loading in layers, and displaying them, but not the style editor yet, although that looks straight forward enough.
However I'm at a loss as to how the variable data should be stored.
For example, while experimenting added a base map layer to QGIS and added a new layer with a  few points on it. I was able to add some parameters to those points, such as an Id and a state boolean.
In QGIS I could add rules then, like if state is true set the point red, otherwise set it green and they would work as expected.
However, in my situation, where I have a database full of node Id's and another database full of GPS co-ordinates, I'm not sure how to format that information in order for GeoServer to be able to parse it properly from my PostgreSQL database.
At first I thought that the points would be stored in geometries and the related dynamic data would be stored in another table linked with an ID or in the same table in an new column.
In order to try deduce this further I exported my points from QGIS in GeoJSON format and got this:
{
"type": "FeatureCollection",
"name": "point_test",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "id": 1, "upgrading": 1, "nodeId": 123456 }, "geometry": { "type": "Point", "coordinates": [ 13.09283649916585, 55.808940217262837 ] } },
{ "type": "Feature", "properties": { "id": 2, "upgrading": 0, "nodeId": 654321 }, "geometry": { "type": "Point", "coordinates": [ 13.030371871652321, 55.959245727217265 ] } },
{ "type": "Feature", "properties": { "id": 3, "upgrading": 0, "nodeId": 987654 }, "geometry": { "type": "Point", "coordinates": [ 13.166037234533263, 55.924109374240864 ] } },
{ "type": "Feature", "properties": { "id": 4, "upgrading": 1, "nodeId": 456789 }, "geometry": { "type": "Point", "coordinates": [ 13.167989254143061, 55.914349276191871 ] } }
]
}

Which has the dynamic properties embedded in the JSON itself. However if I export the layer as an Esri shapefile then I get multiple files, and the properties are separated out in a separate file, in what looks a little more like a traditional database dump, along with some other files. That export looks like this on the drive:

With a few of them not being human readable.
I can import these individual layers into GeoServer no problem, but with the large amounts of data I'll be using it's not practical to add each point in QGIS and export them for use in GeoServer. I'm imagining I'll need to build a service that queries the information I need from multiple databases then builds point objects or geometries from those various sources then dumps them into a postgresql database for GeoServer to use.
So how does something like that normally look in a database? Will those dynamic properties be embedded in a json like structure and queried from there? Or will the dynamic data be stored elsewhere with some sort of ID tying it to the geometry/point?
With the shapefiles, looking at the file ersi_point_test.dbf I can see the dynamic properties stored in there, (and when I click on a point on the map served by GeoServer it displays the correct information) but I can't see a relationship between the geometries/points and the dynamic data; it's obviously located in one of the non-human readable files.


texts:
I need a way to find out the unit of an EPSG-Code (e.g: meter) from java (ogr, gdal, osgeo?) but any other approach would be fine as well.
I want to do something like:

if(epsg.unit==meter){
   do this
} else {
   do that
}


texts: Clustering only if name attribute is same using QGIS?

texts:
I have a point layer on a world map with several attributes. Each point marks the position of a berth at a terminal in a port in a certain country. When I zoom out, it's getting too crowed so I managed to show only the port names and not the terminal and berth names any more by using an expression for the label. This gives me less text so it's already much better. But as a result, I'm now seeing a lot of points all showing the same port name in a certain area which is also kind of suboptimal. I would like QGIS to cluster these points but only if the attribute port name is the same. The individual labels need to be hidden in case of clustering. However, if there is only a single port (unique port name) the clustering is obviously not working but I still need the port name to appear.

I was playing around now for a long time with the Point Cluster function under Symbology but cannot get it to work. Especially the dependency of individual point labels and cluster labels is problematic.


texts: HttpRequestHandler in ArcGIS SDK Android 100.X?

texts:
Il have a service that is exposed to external requests with the adress https://proxyUlrPrefix/arcgis/rest/... It works well to access it and I can display my maps on my mobile app.

When I try to generate a geodatabase to get for offline use though, my server sends the official URL for him on where to get the geodatabase. He gives me an URL like http://myServerName:6080/arcgis/rest/... . Because this URL is not accessible, I cannot get my offline map. So when the request is sent from my mobile to the server, I would like to replace the url myServerName with proxyUrlPrefix so the access is granted.

I saw on the .net SDK that the event ArcGISHttpClientHandler.HttpRequestBegin Event can be used. In the Android SDK v.10.2.7, it's also available as mentioned in Setting up HttpRequestConfiguration to ArcGIS Android SDK for security . In the Android SDK 100.7, I only found RequestConfiguration (ArcGIS Runtime SDK for Android 100.7.0)  that I can use to set headers, but not in order to manipulate the request URL.

Has anyone an idea on why this class is not implemented in the Android SDK and how I can modify my URL?

This question is also posted on GeoNet


texts: Cleanly splitting polylines / reordering vertices in QGIS?

texts:
I'm working with some shapefiles I didn't create and very frequently, when I try to split some of the polylines using the Split Features function, I am left with disjoint parts and more than two segments.

I believe this is related to the vertices being "out of order," but how do I clean up the ordering?


texts: Error on pyqgis dot density maps: no attribute bounding box

texts:
I have been using this tutorial to try to make a dot density map using python within QGIS, and my code shamelessly draws from the linked tutorial:

import random
src = "~/maps/SAMS_areas.shp"
tractLyr = QgsVectorLayer(src, "SAMS", "our")
popLyr =  QgsVectorLayer('Point?crs=epsg:3006', "counts" , "memory")
i = tractLyr.fieldNameIndex('counts')
features = tractLyr.getFeatures()
vpr = popLyr.dataProvider()
for feature in features:
  pop = feature.attributes()[i]
  density = pop
  found = 0
  dots = []
  g = feature.geometry()
  minx =  g.boundingBox().xMinimum()
  miny =  g.boundingBox().yMinimum()
  maxx =  g.boundingBox().xMaximum()
  maxy =  g.boundingBox().yMaximum()
  while found < density:
    x = random.uniform(minx,maxx)
    y = random.uniform(miny,maxy)
    pnt = QgsPoint(x,y)
    if g.contains(pnt):
      dots.append(pnt)
      found += 1
  geom = QgsGeometry.fromMultiPoint(dots)
  f = QgsFeature()
  f.setGeometry(geom)
  dotFeatures.append(f)
vpr.addFeatures(dotFeatures)
popLyr.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([popLyr,tractLyr])

I keep getting an error saying that none type object has no attribute bounding box, and I checked this other question where it is suggested that maybe the cause of such error is the format of the variables. I can confirm that my variables for counts are integers (I use a .csvt file and all). This error appears also if I try to use the dialogue box on the GUI. I wonder if this is because my shapefile has features that are represented by multiple polygons. I ended up with this setup because wanted to spread the dots only within the residential areas of each district/feature, what led to some districts being split in several polygons.


texts: How to fix a calculate value tool to work in a large model and skip featureless rasters?

texts:
I'm attempting to run a model that contains some rasters that do not have any features after being processed through my raster calculator. I have a calculate value tool from a previous question that works exactly as it is shown. The issue is that tool doesn't work when placed within the model. I changed the Data Type to Variant and the calculate value tool runs but returns this error:

The point of the calculate tool was to find these "bad" rasters and skip over them. The calculate value tool correctly identifies and skips rasters in a less complicated model.

Here is the calculate value:

Does anyone know how I can make this calculate value tool work in my big model?

I'm more than open to changing the calculate value code block to run a different process if using Getproperties is causing the issue.


texts: Why doesn't copying values from MapInfo info tool boxes work properly?

texts: As above. When I press ctrl + c, it doesn't always copy the value, I usually end up pressing ctrl + c 5 times to ensure it copies. Is there a way to fix this?


texts: ArcMap v10 - Polygon Layer - Keeping one polygon above another

texts:
I have a shapefile with polygons of differing colors.  I am using a white hash pattern to represent partial ownership.

In short - Polygon #1 is Green and Polygon #2 is a White Hash pattern.

Problem: The Hash Pattern is not staying on "top" of the Color Pattern.

What do I need to do to keep the hash pattern on "top"?

Thanks


texts: R: Using Inverse Distance Weighting with the mapply function

texts:
I would like to find what are the optimal power parameters and number of neighbours for Inverse Distance Weighting (IDW) interpolation using the mapply command and the gstat package in R. The optimal parameters are determined through the leave-one-out cross validation. Until now what I have used for the meuse dataset is the following; for 3 to 10 number of neighbours and for 2 to 5 power:

list_of_idw_results = mapply(idw, nmax = 3:10, idp = 2:5,
                         MoreArgs = list(formula = zinc ~ 1,
                                         data = meuse,
                                         newdata = meuse.grid))


texts:
Is there anyway you can allow users to edit the geometry of shapes (e.g. move an existing point) using the CartoDB.js or CartDB SQL API libraries.

I realise you can do this directly in the CartoDB app but I don't want to give the client access to the CartoDB app and would rather them edit the data through a separate web app.

The only thought I have had so far is to convert the features to GeoJSON and use available leaflet functions to move the data and then update the location via the SQL API.


texts:
I have six ASCII files in ArcMap 10.3 that I am trying to convert into one Raster file for the purposes of producing a linear cross-section of the height (Z) values. I have converted all six ASC files to GDB files, and also to LYR files.

However, when I try the Merge tool to create the final conjoined Raster image, I am getting the error message:

  ERROR 000840 The value is not a Table View.

as soon as I enter the first Raster file to be merged.

If this is not possible, the other possibility may be Raster Calculator, but I am unsure how to proceed with a merge here.


texts: ASTER mineral indices with VINR and SWIR data

texts: I need to run a supervised classification on ASTER data in order to map a rock formation, which is composed of chert. There is no shortage of indices which can be used to highlight silica rich formations, see http://www.ga.gov.au/image_cache/GA7833.pdf. However, all of the indices I have found thus far utilize the TIR bands, which have a spatial resolution of 90m, and I need to map at a much finer scale than this permits. The either of the the 30m and 15m resolutions of SWIR and VNIR respectively should suffice. The problem is that I can't find any silica indices which use these bands. I understand that any such indices might be sub optimal in comparison to those which use the TIR bands, but in this case resolution is of uptmost importance. So, the question is; does anyone know of any established chert/silica indices which use the VNIR or SWIR bands?


texts: Python - OGR: Transform coordinates from meter to decimal degrees

texts:
I have a point in a shapefile. The shapefile is in a WGS_1984_Web_Mercator Projected Coordinate System. I want to get the coordinates of that point in decimal degrees with GDAL/OGR in Python.

I am able to get the coordinates in meters with the code below.

I tried to transform the shapefile from a Projected Coordinate System to a Geographic Coordinate System without success. I assume it is not possible to remove a Projected Coordinate System with OGR. Is that right?

What other ways are there within OGR/Python to get the coordinates in decimal degrees?

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')

shp = driver.Open('testpoint.shp', 0)
lyr = shp.GetLayer()

feat = lyr.GetNextFeature()
geom = feat.GetGeometryRef()

x = geom.GetX() #Get X coordinates
y = geom.GetY() #Get Y cooridnates

print x,y


texts: GEE check only one sentinel tile instead of checking all of them when the featurecollection is in many tiles

texts:
I have code in GEE that creates image collection based on geometry, calc for each polygon the mean value of NDVI for each image, creates new featurecollection with the image as title, and then the polygons names + their mean NDVI value.

My problem is that for some reason the data calculated only for one tile,
originally. I have 10 features:

But in the results, where each polygon name gets column, I get only few (arc, csv and find are the names)
:

I have tries to select manually other tiles, and when I did it, I could get the data of other polygons, but I would like to get all this data at once.

Here is link to my code, I shared the geometry, let me know if it doesn't work:
https://code.earthengine.google.com/b6480f7aa62bbf6ef0736066e095bf68

My endgoal: to get the data for all the polygons, even if they are not in the same tile (image). (or in another words: to create featurecollection with more than one tile-images)

Edit: following Daniel's answer, I have tried to change it from set to create dictionary, so I have tried to change this:

// Collect region, image, value triplets.
var triplets = filter.map(function(image) {
  return image.select('NDVI').reduceRegions({
    collection: geometry.select(['Name']),
    reducer: ee.Reducer.mean(),
    scale: 10,
  })
    .map(function(f) {
      return f.set('time', (ee.Date(image.get('system:time_start'))).format("YYYY-MM-dd"));
    });
}).flatten();

to this:

// Collect region, image, value triplets.
var triplets = filter.map(function(image) {
  return image.select('NDVI').reduceRegions({
    collection: geometry.select(['Name']),
    reducer: ee.Reducer.mean(),
    scale: 10,
  })
    .map(function(f) {
      return ee.Feature(null,ee.Dictionary({time:(ee.Date(image.get('system:time_start'))).format("YYYY-MM-dd")}));
    });
}).flatten();

print('triplets', triplets);

But then when I print "triplets" I get many images that have geoemtry "null" and they have only one column in proprty- time, but without the mean NDVI for the polygon,
 and in addition I get this error message:

  FeatureCollection (Error) Error in
  map(ID=20180911T050649_20180911T052326_T43PHT_00000000000000000000):
  Dictionary: Element at position 0 is not a string.

that if I understand correct- means that it doesn't recognize the column "Name" anymore.

Here you can see the the "modified" code with the featurecollection instead of the set. :
https://code.earthengine.google.com/1df1a654d9473d02c071e28f40d6b014

What am I missing?


texts: What to do when raster cell size is larger than zonal vector feature in order to collect statistics?

texts:
I am trying to extract the mean cell values for a raster dataset based on a survey grid polygon feature using ArcINFO 10.1. The cell size of the raster dataset is about 4 times the size of the survey grid.

I am working with the FCLIM climate trend analysis data for equatorial Africa (ftp://www.earlywarning.usgs.gov/fews/pubs/mapping_decadal_variations.pdf)

I want the mean value of raster cells/pixels within each polygon grid cell:

I have around 3000 grid cells, using 'zonal statistics as table' tool I get only around 600 results evenly spread about the survey grid. Obviously not the result I am seeking.

When I use the 'zonal statistics' tool the result I get the following raster:

2 things have gone wrong with this output (that may not be immediately evident in the image): the cell size of the output raster is still about 4X larger than the feature zone, survey grid polygons; and secondly the raster cells are not coincident with the feature zone survey grid polygons and indeed no longer coincident with the original raster dataset from which it was derived!

My question is: How can I calculate the mean value of raster cells/pixels residing within the smaller survey polygon grids? Is it a matter of using a different method of breaking up the raster cells or is there another method I should try?


texts: How to duplicate this QGIS workflow with gdal/GRASS? (reproject and crop a raster)?

texts:
I have a number of rasters (rgb images) that I want to re-project and crop. I have established a workflow for doing this in QGIS, but I'd like to do it using gdal modules, because the output is being imported into a GRASS location anyway and I could more easily automate the workflow if I were using gdal for the whole process.

Here's my QGIS workflow:
I have a QGIS project that with the CRS (EPSG:32142) set to the target CRS and the desired extent bookmarked (xmin=-985631.39, ymin=7464178.33, xmax=2296529.53, ymax=9628782.84).

I import the raster that I want to convert, zoom to my extent, choose save layer as . . .
In the 'save as' dialog box, I choose:

output mode=raw data
extent=map view extent
resolution=layer resolution

I understand how to use gdalwarp to reproject and set the extent of the output like this:

gdalwarp -t_srs EPSG:32142 \
  inputfile.tif \
  -te xmin ymin xmax ymax \
  targetfile.tif

I also understand (thanks to MappaGnosis) that I can set the resolution for the target file using the -tr argument to gdalwarp. I don't know how to get the source file's native resolution programatically as MappaGnosis suggested in my previous post, and what I'm imagining is a messy process, like using awk or Python to parse the output of gdalinfo to create the required -tr argument for gdalwarp.


texts: How can I transform Client to Server Functions in Google Earth Engine

texts:
I'm trying to transform client-side functions to server function in GEE. However, I cannot obtain the desired result for the function named "add to list". The problem is that I want to add new elements from a list to another list one by one. But the result that I'm obtaining is the first list and then the other one.
https://code.earthengine.google.com/01bddb3d2f7c2eb87d30496646c06382
// Add fixed vallues to a list CLIENT SIDE
var addToList_ClientSide = function(list, season){
  var ini_season = season[0]
  var end_season = season[1]
  var years_seq = [];
  for (var i = 0; i < list.length; i++){
    var el = [list[i], ini_season, end_season]
      years_seq.push(el)
  }
  return years_seq
};

// Add fixed vallues to a list SERVER SIDE
var addToList_ServerSide = function(list, season){
  var ini_season = season.get(0)
  var end_season = season.get(1)
  var years_seq = ee.List.sequence(list.get(0), list.get(-1))
  var el = ee.List([ini_season, end_season])
  return years_seq.add(el)
};


texts: What is the difference between nominal, categorical, ordinal, interval and ratio measurements?

texts:
Given the following examples, which are examples of nominal, categorical, ordinal, interval or ratio measurements?

bushels of wheat per county
vegetation type
average maximum daily temperature
parcel street address
parcel ID number
pH measurements of a stream
state rank for average wage
number of voters in a district
student grade in a class
soil type

In particular, I am having trouble figuring out which one is a ratio measurement.