texts: What is the type of EPSG of longitude/latitude (not Point) co-ordinates?

texts:
I want to know that what is the type of projection (900913, 4326 or so on..) of long/lat co-ordinates. I happen to come across a slight confusion while querying planet_osm_lines when I dump (St_DumpPoints) LineString into point and trying to get long/lat from that point. In the following code, I tried to transform into longitude using 900913 and 4326 projection. I realized that 4326 give me correct result while 900913 give the same as point value.

So, my question is, is 4326 default projection value for long/lat? i.e. any long/lat data that I get from, such as smartphone, or other device are in form of 4326?

In addition, while computing (st_distance, st_closestPoint etc.) the lat/long co-ordinates with the database 900913 based co-ordinates do I have to change it into 900913 projection?
such as:

ST_Transform(ST_SetSRID(ST_makePoint(long, lat),4326), 900913);

COde:

osm=# select d.osm_id, d.name, st_x(st_transform(st_setsrid(d.a,900913),900913)), st_y(d.a) from (SELECT osm_id, name, (ST_DumpPoints(st_astext(way))).geom as a FROM planet_osm_line limit 1000) as d;
      osm_id   |       name        |    st_x    |    st_y
    -----------+-------------------+------------+------------

 328454927 |                   | 9310610.35 | 3223602.19
 328454927 |                   |  9310606.6 | 3223579.35
 328454927 |                   | 9310597.74 | 3223554.21
 328454927 |                   | 9310590.81 | 3223541.52
 328454927 |                   | 9310591.28 | 3223530.75
 328454927 |                   | 9310593.43 | 3223519.49
 328454927 |                   | 9310593.43 | 3223495.79

osm=# select d.osm_id, d.name, st_x(st_transform(st_setsrid(d.a,900913),4326)), st_y(d.a) from (SELECT osm_id, name, (ST_DumpPoints(st_astext(way))).geom as a FROM planet_osm_line limit 1000) as d;
  osm_id   |       name        |       st_x       |    st_y
-----------+-------------------+------------------+------------
 328454927 |                   | 83.6386358188641 | 3223602.19
 328454927 |                   | 83.6386021320409 | 3223579.35
 328454927 |                   | 83.6385225413067 | 3223554.21
 328454927 |                   | 83.6384602880575 | 3223541.52
 328454927 |                   | 83.6384645101394 | 3223530.75
 328454927 |                   |  83.638483823918 | 3223519.49
 328454927 |                   |  83.638483823918 | 3223495.79

I would be very helpful to have a clear understanding on these confusion from you guys.


texts: Ignoring overlapping polygons within same feature class when using Intersect tool of ArcGIS Desktop

texts:
I am attempting to intersect two feature classes, one of which has many overlapping polygons. I don't want to calculate the intersection of overlapping features within the same FC, but instead only the intersection of features between the two FCs.

My problem is similar to Intersecting overlapping polygons individually in ArcGIS Desktop, but the solution proposed there will not work in my case (I'm not trying to calculate zonal statistics.

I've thought of two potential solutions but they each have their drawbacks. First was to simply dissolve the output of the intersect based on the FID field of the FC with overlapping polys. However, this create problems with preserving the format of the attribute table which I need for future calculations

My other idea was to export each individual feature to a temporary FC, then do the intersect and add the intersected features to a new output FC. I think, however, that this would be very slow (i'm doing this over 10,000s of features w/ arcpy).

Is there another way to do this?

(EDIT) Example:
The polygons outlined in red are from the feature class w/overlaps.

When I do an intersect w/ the purple polygons, the result has a 'sliver' where the two red polys overlapped (highlighted below). This is what I DO NOT want.

In this example, the only intersection I really care about is the one highlighted below:


texts: Is there projection issue with Natural Resources Canada (NRC) Canada Atlas hydrology data?

texts:
I have been using the Natural Resources Canada (NRC) data for various generalized maps but it is not until now when i overlaid higher resolution features obtained from the provincial government (Ministry of Natural Resources Land Information Ontario) that I notice a shift in the the NRC Hydrology (watercourse, water bodies) or approximately 800m.

What is strange is that all other data from the same source such as transportation, places, etc. match perfectly with when overlaid with the higher resolution provincial data.

The Provincial data is in NAD83 UTM16 and the Atlas Canada data in WGS84.

ArcMap is projecting on the fly which is usually not an issue (other than performance).


texts:
How can one get image (band) width and height on the server side?
On the user side you do it with image.getInfo().bands[0].dimensions, following Extracting information from Google Earth Engine image by row and column position?
However, it is not possible to map such a function, since getInfo() is a user-side function.
The task sounds extremely simple, since you always see band dimensions whenever you print an image, but how to get them as a list or a numbers?


texts: Choosing method for labeling graduated data that is space efficient?

texts:
Could someone recommend the best standard approach to labeling graduated data ranges?

I was thinking of the following but would like to know if this is industry standard:

Less than 0
>0 - 100
>100 - 200
>200 - 300

instead of

Less than 0
0.01 - 100
100.01 - 200
200.01 - 300

Currently programs such as QGIS default to:

...
0 - 100
100 - 200
200 - 300

Apologies if this is a really simple question but I just can't seem to decide what's best and haven't been able to find any material online.


texts: How to create FeatureLayer without Legend in Table Of Contents (TOC) using ArcGIS Engine?

texts:
So what I am trying to accomplish is having a IGroupLayer that contains an IFeatureLayer and in the table of contents (TOC) there would be no + next to the feature layer (i.e. there is no legend). So far what I tried to do is this, which is not working:

ILegendInfo legendInfo = geoFeatureLayer.Renderer as ILegendInfo;
ILegendGroup legendGroup = legendInfo.get_LegendGroup(0);
legendGroup.Visible = false;


texts:
I would like to convert the tree inventory shape files of London (Canada) into a GPX file and if possible into a JSON file. You can download the archive from their Open Data website. It contains the following files.

READ_ME.txt
TREEINV_SPECIES.xml
TREEINV_SPECIES.xsd
TREES.DBF
trees.gpx
trees.sbn
trees.sbx
trees.shp
trees.shp.xml
trees.shx

Today, I stumbled into ogr2ogr and tried the following command.

$ ogr2ogr -f GPX trees.gpx trees.shp

Though, the following error occurs.

ogr2ogr ERROR 1: Latitude is invalid. Valid range is [-90,90].
This warning will not be issued any more

I already found a similar question but cannot figure out what the problem is. Please mind, I am totally new to the topic of converting geographic data. If you know any better tool I am willing to try it out.

EDIT: An initial success story.

According to the answers of Aragon and Brent Edwards I came up with the following commands. I found a webpage stating that data can be ordered in NAD83. So I assume that the dataset I am trying to convert comes as NAD83 as well. However, I cannot understand how I can find out the EPSG resp. SRID myself. With the given information I can lookup  EPSG:26917 in the EPSG Geodetic Parameter Registry - but if I have no clue what the number is - how would I do that?

$ ogr2ogr -f GeoJSON -s_srs EPSG:26917 -t_srs EPSG:4326 trees_coordinates.json trees.shp
$ ogr2ogr -f GeoJSON -s_srs EPSG:26917 TREES_properties.json TREES.DBF

The converted GeoJSON files contain the following data.

// Example rows from trees_coordinates.json
{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -81.233494, 42.994146 ] } }
,
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -81.232346, 42.994477 ] } }
,
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -81.232176, 42.994527 ] } }
,
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -81.245137, 42.971636 ] } }
,
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -81.189495, 42.996498 ] } }
,
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -81.189651, 42.996442 ] } }

]
}

And ..

// Example rows from TREES_properties.json
{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "TREEID": 19.000000, "DIAMETER": 53.000000, "SPECIES": "MASU" }, "geometry": null }
,
{ "type": "Feature", "properties": { "TREEID": 23.000000, "DIAMETER": 30.000000, "SPECIES": "MANO" }, "geometry": null }
,
{ "type": "Feature", "properties": { "TREEID": 24.000000, "DIAMETER": 53.000000, "SPECIES": "MASI" }, "geometry": null }
,
{ "type": "Feature", "properties": { "TREEID": 104920.000000, "DIAMETER": 5.000000, "SPECIES": "OAWH" }, "geometry": null }
,
{ "type": "Feature", "properties": { "TREEID": 60358.000000, "DIAMETER": 34.000000, "SPECIES": "MAMA" }, "geometry": null }
,
{ "type": "Feature", "properties": { "TREEID": 60355.000000, "DIAMETER": 51.000000, "SPECIES": "MASU" }, "geometry": null }

]
}

Questions:

Are my commands correct?
How can I match the rows of both files? Which property set belongs to which GPS position? Is their a way to combine both in one file? I saw another format that includes all information in one file.
How can I find out the EPSG/SRID and SRS in general if I have nothing more than those files listed in the beginning of my question?


texts:
I'm working with a considerable amount of points and I would like to be able to display only part of them since they overlap a lot and are cumbersome to display.

Most of them form a trajectory so they don't need to be all displayed, depending on the zoom level it is useful to display let say 1 point every 10.

I looked in to the rendering options, I saw you can light or darken the points but I didn't find something useful.I also saw options to display layers only starting form some zoom level (Scale Dependent Rendering as explained in the doc here, but this needs to create different layers of resolution which is not ideal.

Is there a way to solve this problem in QGIS without needing to filter or process the data points in the database before importing them to QGIS?


texts: How to convert poly-line contour feature class into heat map displaying heights

texts:
I'm looking for a way to be able to convert contour data (as lines) into a heat map that would display the data warmest at the highest points and coolest at the lowest.

I have access to QGIS, ArcMap (only the Basic license though), DNRGPS and can use freeware.


texts:
I am having a hard time loading a vectorlayer from the console. Not from the Pythonconsole within QGIS, but from the console on my computer. I later want to write a standalone script!

That is what I am doing

import os
import sys
# where is the module
sys.path.append("/Applications/QGIS.app/Contents/Resources/python/")
from qgis.analysis import *
from qgis.core import *
from PyQt4.QtCore import *

qpp = QgsApplication([], True)
QgsApplication.setPrefixPath(r"/Application/Qgis.app/Contents/PlugIns", True)
QgsApplication.initQgis()

Then I simply to this

vLayer = QgsVectorLayer("/PATH/TO/SHP/test.shp", "test", "ogr")

When I type vLayer I get

  qgis._core.QgsVectorLayer object at 0x110735c30

When I try this vLayer.isValid() I get

  False

Any ideas why?


texts: Using ArcPy to export multiple layers to .pdfs?

texts:
I have found a solution to the issue below.  The problem had to do with coding GetParameterAsText without actually using it as an input in a tool.  That was the reason the raster layers wouldn't turn on.

I am fairly new to Python.  Here is my question:

I have individual range maps for several hundred invertebrate species, each of which has the same extent.  I would like to export these maps to .pdf using Python in ArcGIS 10.2.  Each map also needs to have a "Places" and "Rivers" shapefile (these two shapefiles are the same for all of the species).

I have tried modifying the script from Exporting each layer in map to separate image using ArcPy?, but I can't seem to get the range map rasters to show up.

The following script needs work (and I must give credit to Roy, who answered the other post I mentioned), but this is what I have so far:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]

allLayers = arcpy.GetParameterAsText(0)
lyrList = allLayers.split(";")

PDFPath = arcpy.GetParameterAsText(1)

for lyr in arcpy.mapping.ListLayers(mxd, '', df):
    for layer in lyrList:
        if lyr.name == "Places":
            lyr.visible = True
        if lyr.name == "Rivers":
            lyr.visible = True
        if lyr.name == layer:
            lyr.visible = True
            arcpy.mapping.ExportToPDF(mxd, os.path.join(r"C:\Project_7\Newer_data\Inverts\GO", PDFPath, lyr.name + ".pdf"), resolution=150) #EDIT: This line is working now.
            lyr.visible = False
    arcpy.RefreshActiveView()
    del mxd

I don't think I have added the loops in the correct place.  I need the "Places" and "Rivers" layers displayed on the map every time (i.e., for each range map).

I am unclear how to modify the following to fit my data:

arcpy.mapping.ExportToPNG(mxd, PDFPath+"\\" + lyr.name + ".pdf")

Specifically, what do I replace the "\" with?  Do both backslashes need to be there? (i.e., PDFPath+"C:\Project_7\Newer_data\Inverts\GO\)

Also, I have tried just running the following code without the addition of the two shapefiles, but the layers don't seem to turn on.

for lyr in arcpy.mapping.ListLayers(mxd, '', df):
    for layer in lyrList:
        if lyr.name == layer:
            lyr.visible = True
arcpy.RefreshActiveView()


texts: Getting Only Distinct Values Using Query in ArcGIS iOS

texts:
Can any one Please tell me how can i Get only Distinct Values While Using Query and QueryTask in ArcGIS iOS SDK

This the Code i Write for Getting all the Names.

-(void)viewDidLoad
 {
   NSURL* url = [NSURL URLWithString: @"My Server Url"];
   self.queryTask = [[AGSQueryTask alloc] initWithURL: url credential:cred];
   self.queryTask.delegate = self;

   self.query = [AGSQuery query];
   self.query.where = @"OBJECTID >= 1";
   self.query.outFields = [NSArray arrayWithObjects:@"TYPE_EN", nil];
   [self.queryTask executeWithQuery:self.query] ;
 }

- (void)queryTask:(AGSQueryTask *)queryTask operation:(NSOperation *)op didExecuteWithFeatureSetResult:(AGSFeatureSet *)featureSet
{
  self.featureSet = featureSet;
  NSLog(@"%@",self.featureSet.features);
  for (int i=0; i<[self.featureSet.features count]; i++)
  {
    AGSGraphic *feature = [self.featureSet.features objectAtIndex:i];
    NSString *str = [feature attributeAsStringForKey:@"TYPE_EN"];
    NSLog(@"%@",str);
  }
}

- (void) queryTask:(AGSQueryTask*)queryTask operation:(NSOperation*)op didFailRelationshipQueryWithError: (NSError*)    error
{
  NSLog(@"Error: %@",error);
}

Now i am getting all TYPE_EN in 1000 fields Now i want to get only Distinct Values using Query How can i achieve this?


texts: Write CSV table for join self in ArcGIS for Desktop?

texts:
I wonder how I would write a table for a join which I want to perform.
Right now I wrote a .txt File, tab-seperated. The headerline only contains letters.
Here is the problematic table: http://www.filedropper.com/lookuptable
I am working with the shapes provided here: http://www.eea.europa.eu/data-and-maps/data/clc-2006-vector-data-version

Thats how my file looks like:

ID  CLCCode CLCKlasse
1   111 test1
2   112 test2
3   121 test3
4   122 test4
5   123 test5
6   124 test6
...

Rightclick on the .shp. I choose my .txt file and the two needed fields to perform the join (my key and value pair). My key would be numbers in the attribute table of the shapefile which equals the numbers I wrote in my .txt file.
But, ArcGIS says that there are no matching records.

I also tried it with a , as a seperator. I can choose in the dialog join the textfile, but, the field where it should be stays empty/blanc (it basically gets not selected).

Also, there is something I do not understand. When I open my .csv file in ArcGIS I see the "CLCCODE" column.But, I cannot choose "CLCCODE" in the Join Dialog, although it is a column of the attribute table.


texts: Joining and Editing CSV files in QGIS 3.14

texts:
Three questions, I think are very closely related, if not symptoms of the same issue.
I have added a csv file (no geometry) to my project that already contains shp files.

On import of the csv file, field types are auto-detected. Two of my fields are detected as type int, whereas they should be type double. How can I change the type of the foreign key in the csv file.

The csv file appears uneditable. (It does have a unique field). I have tried toggling on edit mode. I have opened the attribute table, but I can't edit or add data, or change field types.

I have set up a relation linking the primary key (from shp file) and foreign key (in csv file). I have then used the relation widget in the form for the shp file. The primary key is type double. The foreign key is type int. Despite this mismatch, the relationship seems to work when you open the attribute form. But as above, I can't edit or add records.

As an attempted solution, I have tried using "Join attributes by field value" in the processing toolbox, but this fails to join any fields - I suspect because of the int:double mismatch on the join.
Does anyone know how to make the csv file editable and fix the int:double mismatch of the join?


texts: How do we georeference a toposheet using GCS or PCS?

texts:
I am new to the GIS domain. I need to georeference a toposheet, Datum - WGS84, 1:50,000 scale map. I have done the georeferencing using WGS for GCP's and then projected to WGS_94_UTM Zone 43N, using Project Tool.
Another way I followed is, using an online GCS to PCS converter and choosing those values for Ground Control Points. I would like to know which method is correct.


texts: What are the geotiff limitations of geotools?

texts:
I'm trying to process globcover, which is a .tif file of 374mb. Running the following code:

File tiffFile = new File("C:\\Users\\CAA\\Desktop\\test.tif");
GeoTiffReader reader = new GeoTiffReader(tiffFile);
grid =reader.read(null);
gridData = grid.getRenderedImage().getData(); //row 32

I get the following error:

java.lang.IllegalArgumentException: Requested region cannot be represented by a single Raster.
at javax.media.jai.PlanarImage.getData(PlanarImage.java:2163)
at javax.media.jai.PlanarImage.getData(PlanarImage.java:2016)
at javax.media.jai.RenderedOp.getData(RenderedOp.java:2266)
at de.myapplication.de.myapplicatoin.GeoTiff.initTif(GeoTiff.java:32)

After playing around a bit with QGis and it's Clipper tool i noticed that if i extract a smaller area, the reader works as intended. If i extract huger areas (i.e. entire europe or the entire world) i get the error again. Maybe it's just correlation and not causation, but this leads me to believe that there is somewhere a limit inside the reader. But what is the limit?

As requested the output of gdalinfo for the original tif file.

Size is 129600, 55800
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 = (-180.001388888888900,90.001388888888883)
Pixel Size = (0.002777777777778,-0.002777777777778)
Metadata:
  AREA_OR_POINT=Area
  Authors=Sophie Bontemps <sophie.bontemps@uclouvain.be>,Eric Van Bogaert <eric.vanbogaert@uclouvain.be>,Pierre Defourny <pierre.defourny@uclouvain.be>
  Copyright=Copyright ┬®UCL Geomatics, BELGIUM 1999-2010
  Generate by=gdal_mean
  process begin time=2010-12-23T09:49:37
  process files=CL5_GLOBCOVER-L5_CLASSIF_2009_V2.3.20101220.tif
  process finish time=2010-12-23T09:57:38
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0013889,  90.0013889) (180d 0' 5.00"W, 90d 0' 5.00"N)
Lower Left  (-180.0013889, -64.9986111) (180d 0' 5.00"W, 64d59'55.00"S)
Upper Right ( 179.9986111,  90.0013889) (179d59'55.00"E, 90d 0' 5.00"N)
Lower Right ( 179.9986111, -64.9986111) (179d59'55.00"E, 64d59'55.00"S)
Center      (  -0.0013889,  12.5013889) (  0d 0' 5.00"W, 12d30' 5.00"N)
Band 1 Block=129600x1 Type=Byte, ColorInterp=Gray
  Min=11.000 Max=230.000
  Minimum=11.000, Maximum=230.000, Mean=183.164, StdDev=55.663
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=230
    STATISTICS_MEAN=183.16435194952
    STATISTICS_MINIMUM=11
    STATISTICS_STDDEV=55.662605157347


texts: Weird differences between basemap rotate_vector and Cartopy transform_vectors

texts:
I define a Lambert Conical Conformal in both Cartopy and Basemap
    import numpy as np
    from mpl_toolkits.basemap import Basemap
    #Cartopy:
    lam2 = ccrs.LambertConformal(central_longitude=-2.5, central_latitude=40.0, standard_parallels=(40.0,40.0))
    latlon = ccrs.PlateCarree()
    #Basemap:
    m = Basemap(width=12000000,height=9000000,projection='lcc',
           resolution='c',lat_1=40.0,lat_2=40.0,lat_0=40.0,lon_0=-2.5)

OK, so far, I think it is correct. Now, let us rotate vectors, in Basemap:
m.rotate_vector(np.array([-2.6257]), np.array([-4.944]), np.array([349.835]), np.array([34.475]))
Out: (array([-3.04061787]), array([-4.70022121]))

But if I do the same operation in Cartopy:
    lam2.transform_vectors(latlon, np.array([349.835]), np.array([34.475]), np.array([-2.6257]), np.array([-4.944]))
    Out: (array([-2.68554888]), array([-4.91174752]))

A result that is significatively different to the result from Basemap (especially in the u component of the vector). Am I doing something wrong? Is this normal?


texts: Aggregation of multiple POSTGIS records in Geoserver

texts:
I have a WMS layer which geoserver fetches from a postgis layer. For a given query to this layer I might have two records for a given geometry. The record contains geoshape and value where value is used to create lengend graphics. How does geoserver aggregate these two records?

An example of 2 records for the same point:

|geometry                          | gridcellname | value |
|----------------------------------|--------------|-------|
|POINT(3.37531652151541 80.7506792582348)|"Y71X20"|0      |
|POINT(3.37531652151541 80.7506792582348)|"Y71X20"|137.756698608398|

Does geoserver
1. take the max value or
2. take the average value?
Does the aggregation happen differently if you have a POINT() geometry verces a MULTIPOLYGON() geometry?


texts: Is there any way to have all the required files of JavaScript API? Because I have a local application and the user won't be able to have the resources from serverapi.arcgisonline.com


texts:
I have a polygon layer (>600000 features) in a Spatialite database and I want to overcome to run out of RAM by selecting only the features I need for subsequent workflow.

I can connect to my Spatialite database using RSQLite with the following code:

dbfile="Path/To/field.sqlite"
sqlite=dbDriver("SQLite")
con=dbConnect(sqlite,dbfile, loadable.extensions=TRUE )

I can also query my table and get an tabular output:

above1600<-dbGetQuery(con,"SELECT * from height_aspect where height >= 1600")

But I do not understand how to get a SpatialPolygon layer.

Do I have to convert the output or do I need to initialise the database as a spatial one?

I researched the web for a while but could not find anything that helps me.


texts: Is schema definition compulsory for WFS-T

texts:
Most tutorials on WFS transaction uses a schema and geom attributes in WFS request like this

var wfs_layer = new OpenLayers.Layer.Vector("Editable Features", {
strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
//projection: new OpenLayers.Projection("EPSG:26910"),
protocol: new OpenLayers.Protocol.WFS({
version: "1.1.0",
url: "http://localhost/geoserver/wfs",
featureNS :  "http://www.opengeospatial.net/cite",
featureType: "wfst_test",
**geometryName: "the_geom",**
**schema: "http://localhost/geoserver/wfs/DescribeFeatureType?version=1.1.0&;typename=cite:wfst_test"**
})
});

can we omit this or is it compulsory to include them?


texts: Displaying pointcloud symbolized with RGB values in ArcMap?

texts:
I know that I can display a LAS dataset in ArcScene with RGB Symbology, but is there a way of doing this in the 3D View window in ArcMap? The LAS Dataset toolbar only has symbology for elevation, class and return. Why is intensity and RGB not available to select like in ArcScene?

I am using ArcMap (version 10.6.1) which includes a 3D MapViewer with the LAS Toolbar.


texts: Selecting intersected lines from two common attributes in very big data set?

texts:
I have a layer with lines representing cables. The data set is huge, over 100 000 polylines. Some of these cables lies within the same shaft, some of them lies within the same shaft as an other cable for a bit and some of them lies within their own shaft.

Each individual line has the attributes of a unique ID, the year it was build, the voltage and the lenght of the line.

What I need to get to know is where these lines intersect with eachother. However, I also need to present cables dug down in the ground at the same time and the lenght of the shaft were they go together.

What I've done so far is to intersect the linelayer and by doing that I get a very very long list of were the lines intersect. My struggle is that if I want to keep all the attributes from the input layer, my output layer is even more fragmented. This is due to if on long cable intersect with three shorter ones, the output consist of six pieces.

What I realised is that I can sort the attribute table in the output layer after lenght and visually se the "shaft groups" since the output consisted of so many decimal no shaft was the same lenght as another. However, the dataset is way to big to sort it manually. So that is pretty much were I'm at right now, with a huge list of cable-intersects and the need to select all cables that are in the same shaft, and the lenght of the shaft, and from the same year.


texts: Vectorising lines appear as a mixture of bubbles and lines from raster map

texts: I have these road raster maps which I have reclassified into binary images to try and vectorise. However, when I click on vectorise and generate feature, it comes out as a mixture of bubbles and lines (see screenshots). I tried to play around with the vectorsing settings but without much success. I would like to just get a line following the roads.


texts:
I'm struggling to find a quick way to identify and label an area by the maximum coverage of an overlapping polygon. Has anyone done a process like this lately?
I'm currently playing with the Pivot table, but the machine is still processing, and I'm not sure if the Pivot Table is the answer, yet.
I'm using ArcGIS 10.1


texts: Raster with wrong extents when using Polygon to Raster

texts: I am trying to use the "Polygon to Raster" tool on a building extents polygon file  to create a raster with the exact extents and cell size as an existing flood depth raster. Even though I have added my desired raster extent in the snap raster option in the Environment settings and defined the cell size, the result does not have the same extents. Furthermore, I made sure that they have the same projection (WGS_1984_UTM). I really need them to have the same extents and cell size as I want to do some calculations between them (e.g. calculate the flood depth at the buildings location). Any idea why this is happening?


texts: How To Display Service Flow Direction Graphics on User Side( JS API or AGO)

texts:
Using ArcGIS Server 10.4.1 I am publishing a Service which has Flow Directions on it like

but when consuming in ArcGIS online or even on ArcGIS for JavaScript API Application I am getting  lines without Arrow


texts:
I'm building simple views against an ArcSDE Oracle RAC 10g. My SHAPE columns are of SDO_GEOMETRY datatype. ArcSDE 9.3.1.

I have several point tables, with the following outline:

ID (Guid)
Name (string)
Shape (SHAPE)

I've built a simple view liek this in Oracle!

CREATE OR REPLACE VIEW VW_POINTS AS

    SELECT ID, Name, Shape from TABLE_A

        UNION ALL

    SELECT ID, Name, Shape from TABLE_B

        -- etc

The view in Oracle is created just fine. I've registered this view with the command line:

sdetable -o create_view -T POINTS -t VW_POINTS -c Id, Name, SHAPE etc

The ArcSDE view is created again, just fine. I can use ArcObjects to query it, etc. I need now to create a Layer. I've tried the command line tools, sdelayer -o register, but I had no success. It gives me a message, saying that views are not suported.

I will need to do this several times, sometimes using complex syntax for building the view.

What is the problem? How can visualize these spatial views as SDE layers (I don't need to edit them, just view them spatially inside ArcMap)

Thanks for the input!

EDIT

This is a screenshot in ArcMap, for your deligth :P

ArcMap Screenshot for your happiness:

As requested, here comes the describe_long

// create_view
C:\Users\H>sdetable -o create_view -s server -i port -u user -p pass -T
view_airspace_points -t vw_airspace_points_oracle -c "*"

at this point I altered the view on oracle, including a cast as numeric(38,0) before registration.

// register layer
C:\Users\H>sdelayer -o register -l view_airspace_points,shape -C OBJECTID,USER -
g RTREE -s server -i port -u user -p pass -e p -P HIGH

ArcSDE 9.3.1  for Oracle10g Build 2784 Tue Oct 27 10:51:14  2009

Layer    Administration Utility

Layer Description ....:

Table Owner ..........: SIGMAGIS2

Table Name ...........: VIEW_AIRSPACE_POINTS

Spatial Column .......: SHAPE

Layer Id .............: 1398

SRID .................: 79

Minimum Shape Id .....: 1

Offset ...............:

falsex:       -180.000000

falsey:        -90.000000

System Units .........:     994200.000000

Z Offset..............:          0.000000

Z Units ..............:          1.000000

Measure Offset .......:

Measure Units ........:

XY Cluster Tolerance .:          2.0

Spatial Index ........:

parameter:    SPIDX_RTREE

exist:        Yes

array form:   -2,0,0

Layer Envelope .......:

minx:      -180.00000,        miny:       -90.00000

maxx:       180.00000,        maxy:        90.00000

Entities .............: p

Layer Type ...........: In-Line Spatial Type

Creation Date ........: 06/28/11 09:28:45

I/O Mode .............: NORMAL

Autolocking ..........: Enabled

Precision.............: High

User Privileges ......: SELECT, UPDATE, INSERT, DELETE

Coordinate System ....: GEOGCS["Longitude / Latitude [WGS 84]",DATUM["WGS 84",SPHEROID["WGS 84",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Decimal

Degree",0.0174532925199433]]

Layer Configuration ..: DEFAULTS


texts: How to trigger a hover event for map graphics in ArcGIS Javascript API?

texts: I want to trigger an event when the user hovers over a geometry that they drew. I know how to set up hover events for DOM objects, but not for graphic elements within the map window.


texts: How to build tiles and display them without georeferencing?

texts:
When I work with georeferenced image (GeoTiff) I can use gdal2tiles.py to build tiles and generate googlemaps.html, leaflet.html, openlayers.html, etc. which I can upload to server and share with other users.

But what if I want to just created tiles of big tiff image (not geotiff image for example image of moon) and share it without georeferencing to any map? How it can be done?


texts: Can't use simple GRASS modules in python (python script not found)

texts:
I am using grass 7.8.4. I'm trying to run some simple commands using the Grass Python Scripting Library in Python. For modules for which the source code is written in C++ or any other language than Python, it works without problem but for the modules written in Python it does not. Here is an example of a simple command I tried to run:
import grass.script as gscript
gscript.run_command("v.db.dropcolumn", flags="help")

It raises the error:
grass.exceptions.CalledModuleError: Module run None v.db.dropcolumn --help ended with error.
Process ended with non-zero return code 9009

I already checked and I know that return code 9009 means it does not find the python script. I solved the error by just adding ".py" in the command:
gscript.run_command("v.db.dropcolumn.py", flags="help")

I would like the command to work without the ".py" because I'm writing a python script that makes use of some addons. For the addons to work in my python script I sometimes need to modify them because there are running commands like "v.db.dropcolumn" without the ".py", so I have to add it. But then, the addons do not work anymore when I use them directly in the grass GUI. Does anybody know how to solve this?


texts: How to register a broadcast on QGIS to start my plugin on start after loading plugins

texts:
I want to perform a default task when user opens QGIS but here I am building it as a python plugin so I want a signal or broadcast to call my plugin.

I didn't build QGIS from source
my plugin should be called when QGIS finish loading plugins on start

can anyone help with the documentation I can refer about the broadcasts and calls i can register to.

how should I approach my problem?
I have found a help on startup.py but it seems to help run python not run plugins.but tried using it and keep bugging errors like

Couldn't load qgis.user.

Python support will be disabled.

Now my plugin is python/plugin/Loadmap folder in side it there is a file LoadMap.py and I want to run the run method inside it.


texts: Halftone effect gridded points inside polygon

texts:
I am trying to create smooth looking polygons by gridded points in QGIS (2.18.16), with a fading effect towards the boundary of the polygons (halftone effect).

By using the 'fillGrid' function from Rob J, I managed to create gridded point pattern in Geometry Generator.

Now I would like to add a size decrease of the point marker towards the boundary of each polygon starting from a maximum size (e.g. 2 mm) to a minimum size (0 mm). What is the best way to do this? Maybe there is a faster way by  using the standard 'point pattern' option as symbol layer type?

Example image demonstrates the desired effect inside each polygon.


texts: Can I have ogr ExportToJson export with NO z coordinate?

texts:
I have code exporting a polygon to a json file, I would like to do this without the z coordinate since I don't need it and editing the database that was created from geodjango migrates looks painful.

So some way to tell this:
poly.ExportToJson()  to not include a z coordinate?

snippet that exports the json from a read in geotiff file:

for gcp in gcps:
    ring.AddPoint(gcp.GCPX, gcp.GCPY)

ring.AddPoint(gcps[0].GCPX, gcps[0].GCPY) #close off the polygon
    # Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
md = nitfData.GetMetadata()
print(filename)
dic = {"sensor" : md['NITF_ISORCE'], "corner_coords" : md['NITF_IGEOLO'], "target_name" : md['NITF_TGTID'], "collection_dt" : md['NITF_IDATIM'], "f_size":os.path.getsize(path), "polarization" : md['NITF_FTITLE'], "mpoly" : poly.ExportToJson() }

jsonPaths.append(json.dumps(dic) + ',')


texts: Stopping ArcMap from automatically changing layer source path after upgrade?

texts: I have just upgraded to ArcMap 10.5. I noticed a change that I do not like or currently want. When I copy an MXD to a new folder and/or move layers from one mxd to another it is automatically changing the layer paths to the home folder. This is causing a disaster in all my maps. I cannot find the option to stop it from doing this.


texts: Get a list of all supported EPSG code from OGR

texts:
I need get a a list of all supported EPSG code from OGR/GDAL. I need this becaouse i'm developing a C++ application map, and I need this information to put it in a user form.

I read the OGR documentation but I couldn`t find any method to do this. Some one know how is possible get this information in OGR?


texts: Add featureLayer as operationalLayer from JSON in CMV

texts:
Is it possible to define operationalLayer as feature (or anything else?) dynamically from json?

We have API on our server which generates JSON as a query result from database (acting as a middleware), I can modify JSON through PHP and I would like to add this to cmv operationalLayers so I can toggle the layer on/off through layerControl. Then I need info popup on this layer and custom symbology based on categories.

I don't know how should I approach this - create KML? But the server is not accessible from internet, so it didn't work out. Create CSV and use csvLayer? Well, it works somehow, but I don't know how to style the layer based on attribute unique values. It looks like featureLayer is the best approach. According to this example  it works on JS API, but I don't know how to implement it in CSV.


texts: Equivalent of "Dissolve with stat" in QGIS graphical modeler

texts:
I'm trying to sum multiple variables by categories (id of my polygons) within a QGIS model but I haven't found a proper way to do that...

If I hadn't to use the graphical modeler I would obviously use "Dissolve with stats" but this tool cannot be added within a model...

Is there another way to do such thing without using the python console ?

Here is a picture of the attribute table of my shapefile:


texts: Save QGIS style created with espression string builder

texts:
I created a rule-based renderer with espression string builder, where I also use the function editor in Python. When I share the style (qml) file with others, I realized that the Python function I wrote is not saved in the style file.

Is it possible to save the Python script inside the qml style file? And another question is, where are the Python script saved? I checked the usual pathes where QGIS stores script, but was not able to find it.


texts: QGIS Processing Modeller reusing output as input for next step

texts:
I would like to run the Sample Points algorithm over 20+ rasters, reusing the sampled result as the input for the next iteration of Sample Points.

At the moment I am achieving the wanted result by using the Batch Interface and calculating the output of Sampled Points by using an expression like:

'Sampled_' || @row_number || '.gpkg'

And doing similar thing for the input of Sample Points (apart from first row where original sample point file is referenced):

'Sampled_' || to_string( @row_number -1 ) || '.gpkg'

Is there a better way of doing the same workflow without creating 20+ temporary .gpkg files?


texts: Saving NDVI TIFF file as txt file in ArcGIS Desktop?

texts:
I have to use elevation data (dem.tif) and NDVI data (ndvi.tif) as covariates for regression kriging and I want both these files as txt files.

How do I do this?


texts: Is there a LINQ-provider for GeoDatabase

texts:
I want to query my ArcSDE geodatabase using LINQ but I don´t know if there is an ESRI-provider that may also read geometries and which is able to handle archiving as well. Thus I´d like to read data from a FeatureClass for example and apply a whereclause:

IRelationalOperator relOp = (IRelationalOperator) myFeature.Shape;
var features = featureClass
    .Where(x => relOp.Within(x.Shape) && x.MyProperty > 1)
    .OrderBy(x => x.MyProperty)

Appaerently while querying simple structures such as a numeric field might be quite easy I suppose geometric queries expect much more work.

Does this even exist or has anyone ever done something similar? Will it exist on ArcGIS Pro (of course not using arcobjects)?


texts: R shapefile transform longitude latitude coordinates

texts:
I want to find out the coordinates (latitude and longitude) of Texas precincts. I download the shapefile from  here.

I am using R to read the file and extract the coordinates:

dsn.tx <- "~/.../Texas.shp"
shp <- readOGR(dsn=dsn.tx, stringsAsFactor = F)
coords <- coordinates(shp)

But the resulting coordinates are in an unusual format:

    X         Y

 1411246  864457.9
 1266383 1231047.3
 1213326  905663.0
 1307353 1213580.7
 1168558  814300.4

I am not sure why this happens. In other (not sure if related) posts some suggested it may be related to the proj4.

How do I transform the resultant coordiantes into "standard" longitude latitude coordinates?


texts: Adjust page size to map size automatically in QGIS print composer

texts:
I'm looking for a stable way to create a print composer template which automatically adjust its size to the area of interest (defined by a layer "perimeter" containing one polygon feature) at fixed scale 1:25'000.

I tried using the Atlas tools and set it to fixed scale. This works, but it's not stable and difficult to handle with different projects (different perimeter layers).

It most often does not render the map when Atlas is set to fixed scale in "controlled by Atlas", which means changing to "margin around feature". Further, it always kicks my preset scale and uses "nan", and I haven't found the logic to recreate it.

Any suggestions?


texts: Why is i.histo.match missing from QGIS 3.x?

texts:
I'm trying to normalise and then mosaic two ASTER scenes (multiple bands in each). The GRASS tool i.histo.match has been suggested in a few places, but it isn't listed in the processing toolbox for QGIS 3.x

https://grass.osgeo.org/grass74/manuals/addons/i.histo.match.html

Can anyone advise how to have this script added to the GRASS tools of QGIS?


texts: QGIS 2.2 crashes when saving digitized data in ArcSDE 10.1

texts:
I am using QGIS 2.2 with ArcSDE 10.1 for Oracle 11gR2
I have created a polygon feature class with SDO_geometry, when trying to save the digitized data QGIS crashes. Please help, I want to edit the ArcSDE data in QGIS.


texts: ArcMap document "Accessed" timestamp incorect display in Windows Explorer

texts:
I created an ArcMap document on February 22, 2016. The timestamp, as seen in Windows Explorer, is 2/22/2016 8:58 AM.

I have updated and saved that document many times since, but the timestamp doesn’t change. When I look at the document properties in Windows Explorer, I see that Created (2/22) and Modified (today) are correct, but Accessed is the same as Created, instead of Modified.

This peculiarity affects this MXD only. Other MXDs behave as expected (i.e., "Accessed" matches "Modified"). Any ideas?

I am running ArcGIS Advanced 10.3.1 on Windows 10. The document is stored locally on my C: drive.

I reached out to Esri tech support first, who gave me a superficial, wrong answer. So here I am.


texts:
Reading the ArcGIS help: Supported raster dataset file formats, I read that ASCII Grid file types are addressed with the singular file extension .asc for both 16 bit signed-integer (for discrete data) and 32 bit floating-point files (for continuous data).
They say in ESRI Grid format section:

A grid is a raster data storage format native to ESRI. There are two types of grids: integer and floating point.

and then, it describes what both formats are.
My confusion comes in the next help section ESRI ASCII raster format. There, they only mention the floating point format. So, I would like to know:

Is the 16-bit signed integer format for discrete data also considered a raster format, besides grid?
In 'ESRI Grid format - Grid data storage' section, they cite several other extensions such as .bnd, .hrd, etc.? Is the 16-bit signed-integer .asc-file, the ASCII version of an Esri discrete data grid format?

The motivation for this question is to describe/understand what is a file with extension .asc.


texts: How to transfer UTM coordinates to Euclidean coordinate system?

texts:
I get the coordinates of objects in the UTM format:

"coordinates": [
                [
                    4631430.35,
                    396976.54
                ],
                [
                    4631428.01,
                    396972.19
                ],
                [
                    4631430.24,
                    396970.98
                ],
                [
                    4631432.59,
                    396975.33
                ]
            ]

Also, I have the UTM reference => 32N.

How can I transform these coordinates to the Euclidean coordinate system? Which additional data may I need to do it?
I am looking for open-source Python libraries that can do such transformation.


texts: How to add a bounding box to OpenLayers 3 map

texts:
I have built a simple map using OpenLayers 3, have a look at the code below:

JavaScript code:

var map;

function init(){
    map = new ol.Map({
        target:'map',
        renderer:'canvas',
        view: new ol.View({
            projection: 'EPSG:900913',
            center:[-8015003.33712,4160979.44405],
            zoom:5
        })
    });

    var newLayer = new ol.layer.Tile({
    source: new ol.source.OSM()
    });

    map.addLayer(newLayer);

    /*
        get and display extent
        http://gis.stackexchange.com/questions/168590/how-to-get-extent-in-openlayers-3
    */

    var size = map.getView().calculateExtent(map.getSize());

    console.log(size);

    console.log(map.getSize());

    /* ========================== */

}

/* function to remove layer */

function removeTopLayer(){
    var layers = map.getLayers();
    layers.pop();
}

/* ========================= */

/* swap layers function */

function swapTopLayer(){
    var layers = map.getLayers();
    var topLayer = layers.removeAt(2);
    layers.insertAt(1, topLayer);
}

/* =================== */

HTML code:

        <button style='position:absolute; left: 10px; top:30%' onclick='removeTopLayer()'>Remove Top Layer</button>
        <button style='position:absolute; left: 10px; top:40%' onclick='swapTopLayer()'>Swap Top Layers</button>

Functioning link can be found HERE.

Now I wanted to add a bounding box, so that when the map loads, it only loads for a specific area, I searched for examples online, but  they are either of OL2 or are about a user selecting a bound box dynamically, which is not my requirement. My requirement is pretty simple, all I want is the tiles for a specific bounding box to load instead of the map for the entire world to load.

I found the below method in OpenLayers 3.0 doc's boundingExtent, but I am not sure how to use it. :(  So now suppose I want to add a bounding box for THIS area, how would I go about doing it, using the boundingExtent method? Can anybody explain?

I have been scratching my head over this one for quite sometime now, can anybody help me achieve my goal of have only a specific bounding box load?


texts: Change Order of Features in Popup with ESRI ArcGIS Javascript API version 4

texts:
I am using version 4.5 of the ESRI ArcGIS Javascript API with a MapImageLayer.
I'm setting popup templates on each layer.

The problem is that, when showing the Popup following a click in the map, if there are multiple features found, I need to be able to order the features found.

Is there a way, like in the version 3 of the api, to order these features before they are displayed by the Popup?


texts: Is there a service to check whether a gps path belongs to a public road?

texts:
I have a bunch of GPS tracks. I want to exclude the ones which go via public roads. I'm interested in off-road segments only. Is there any public api to check if a given set of coordinates belongs to a road or street? Is there a better way to do it? Perhaps OSM?

Some context - the tracks are in UK. I'm comfortable with programming in C#/F#. I'm thinking of writing a piece of code which will filter the tracks. Computing power is limited to a single awerage PC.


texts: How to find the polygon containing a given coordinate using ArcGIS Geometry Service?

texts:
Say I have an ArcGIS Server GeometryService running, and a map service containing ZIP code polygons (all ArcGIS 10.1).

I want to pass a coordinate, and have it return which ZIP code contains that coordinate. There seem to be many ways to approach "does THIS point fall within THIS polygon?", but how would I find out which polygon among them all contains a point?

In other words, there's this example where, given two points and a polygon, it will return which of the points is contained by the polygon. But what if you only know the point, and the objective is to determine which polygon?


texts: Comparing cities polygon over time using Nightime lights

texts:
I am currently generating polygons of cities in India over time based on night light data. For this, I have been using the clump function in R and then converting them to polygons.
My issue comes when I create the cities for the next year. These new cities do not share any common properties from the cities created the previous year.
How can I compare changes in these cities along time?

For example, as shown in the image, in red are the city boundary in 2000 and in blue the city boundary in 2010. I need to somehow make sure that they share a common property so then I can compare them.


texts: How do you memorize the meaning of longitude and latitude?

texts:
Longitude and latitude are concepts I usually don't need to remember them. However I sometimes suddenly need to learn their meaning for a period of time because I am using them.

Is there an aide-memoire(memory - aid) to finally remember their meaning until the end of my days?

What is the etymology of these two terms?


texts: Stacked label with parentheses and exclusion

texts:
I'm doing some labeling and have my labels stacked using python by using the following: [field1] + 'n\' + [field2]. But I'd like to keep the label stacked and add parentheses to [field2] as well as ">= 1". I tried using the following with little success. So really I don't know how to combine everything correctly to get it to work properly.

def FindLabel ( [field1], [field2] ):
  if long([field2]) >= 1:
    return "[field1] + '\n' + [field2]"
  else:
    return [field1]

Ideally, I want my label to look like this:

Field1
(field2).

Got it to work besides the parentheses. I would still like to add them if possible.

def FindLabel ( [NAME], [Deaths] ):
 if long([Deaths]) >=1:
  return [NAME] + '\n' +  [Deaths]


texts: Change rasters' z value to minus (Below sealeave)

texts: My current raster has positive values for all features but I need them all to be negative as it is for a seafloor basin. A handful of values are negative but most are positive. I have access to ArcGIS Pro.


texts: Overpass API: Get coordinates of postal boundary

texts:
How can I use the overpass-API to recursively get all nodes and their coordinates which belong to the ways of an relation? I need this to draw bounding polygons of all the postal areas inside a city.

I have the following query:

<osm-script>
<area-query ref="3600062649" />
  <recurse type="up" />
  <query type="relation">
     <item />
     <has-kv k="boundary" v="postal_code" />
  </query>
<print />
</osm-script>

This returns all relations with postal boundary inside the queried area and all the ways of these relations, however it does not return the nodes/coordinates of the points which make up these ways:

<relation id="1300371">
   <member type="way" ref="51702433" role="outer"/>
   ...
   <tag k="note" v="04451 Borsdorf"/>
   ...
</relation>
<relation>
....
</relation>

I would like to get a list of postal areas with their nodes, like this:

<relation id="1300371">
   <member type="way" ref="51702433" role="outer">
      <node id="323735839" lat="51.3751873" lon="12.5689796" />
      ....
   </member>
   <tag k="postal_code" v="04451"/>
   ...
</relation>

The tag with the postal code should still be included so I can distinguish all received postal areas.


texts: Categorizing & Legend with two attributes in QGIS

texts: On a Layer, preferably a point, I categorized its style. So after categorizing using attribute1, I would like to show the legends using attribute2. How can I do that?


texts: How to draw a square from 2 coordinates?

texts:
I have a CSV file with the following output

id, x min, y min, x max, y max, value
1 , 0    , 0    , 100  , 100  , 5
2 , 0,   , 100  , 100  , 200  , 10
3 , 100  , 100  , 200  , 200  , 15
4 , 100  , 0    , 200  , 100  , 20

I can plot the co-ordinate points onto my basemap fine. What I'm trying to do is tell QGIS that row 1 of the CSV has the bottom left and top right corners of a square, to draw it, and to assign it a value of 5.

I'm not even sure this is the right way to achieve my goal - which is to display a grid with different values attached to each square, and then produce a heat-map based on those values.

Thanks.


texts:
I have two polylines, the first consists of individual GPS traces and the second is a series of road network poylines.

I wish to match the GPS trace polyline to that of a road network polyline or series of polylines in order to pass the time attribute from the GPS traces to the road network poyline(s). Of course, for a range of accuracy related reasons the GPS trace does not perfectly match the road network.

This may be visually inspected. The blue line shows the GPS trace polyline and the gray represents the road network polylines.

Due to the complexity of road junctions and directionality restrictions, a nearest neighbour match for the range of coordinates within the polyline won't work.

An example GPS trace polyline looks like this:

GPS_trace = [
"polyline_decoded":[
          [51.5201035,-0.0965073],[51.5201,-0.09651],[51.52004,-0.09693],[51.5201406,-0.096996]],
"duration":15]

An example road polyline looks like this:

road_network_link = [
"polyline" : [
[51.5223,-0.0965],[51.521,-0.097],[51.52,-0.0979],[51.5212,-0.0970]]
"id" : "10200A"]

This appears to be the problem a good solution - https://developers.google.com/maps/documentation/roads/snap but I'd prefer to try myself first.

I have a large amount of GPS polylines and I am thus looking for a matching process that scales well. I'm aware this functionality exists in ArcGIS and QGIS.

My first inclination is to perform a shortest path algorithm for the GPS polyline points via the road link network network. Can anyone point me in the right direction?


texts: QGIS 2.0 Problem with loading a Layer Style file that contains a transparent pixel list

texts:
I work with QGIS 2.0.1.

For my raster file I have created a Layer Style with a transparent pixel list.

When I import a rasterfile in a new project and load the layer-style file, the map changed to the correct representation.

But the transparent pixel list is empty.

When i want to apply the style settings, the presentation changed back to the original shape style.

I have looked in the style file and the transparent pixel list are stored.

Code:

<rasterTransparency>
    <singleValuePixelList>
      <pixelListEntry min="0" max="0" percentTransparent="100"/>
      <pixelListEntry min="1" max="1" percentTransparent="40"/>
    </singleValuePixelList>
</rasterTransparency>

Second workflow:

Open the Properties from the raster file
Select the the Transparency menu
Click the + Button (add values button)
Enter the raster values and adjust the transparency
Export the list to .txt file

After this I load the transparency pixel list text file into the new Project additionally to the layerstyle . This works.

Why is the transparent pixel list stored in the layer style file but you must store the pixel list as a extra text file?

What can you do, that you don't have to save 2 files (Layerstyle and transparent pixel list .txt)?


texts: Placing point over pixel with maximum value falling in grid overlay using QGIS

texts:
I have a raster layer with a Shapefile grid with cells measuring 3 x 3 m overlaid on top. All I want to do is find the pixel with the maximum value within each cell, and place a Shapefile point over that pixel such that each cell has a point. See figure below (excuse the poor job).

I've found Q&As that show how to identify the maximum values, but none that transform the query to a point explicitly. The point is that the points sit above the maximal pixel, not the cell centroid.


texts: Trouble opening Sentinel 2 image with pillow

texts:
I am doing data augmentation to create training data but the Sentinel 2 image is not showing. It always gives me an error cannot identify image file.
Does someone has any solution for this or maybe another methods to do data augmentation to create the training dataset.
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
import io
from PIL import Image
# Construct an instance of the ImageDataGenerator class
# Pass the augmentation parameters through the constructor.

datagen = ImageDataGenerator(
        rotation_range=180,
        width_shift_range=0,
        height_shift_range=0,
        shear_range=0,
        zoom_range=0.5,

        horizontal_flip=True,
        vertical_flip = True,
        fill_mode='nearest')    #Also try nearest, constant, reflect, wrap

a = 1
while True:
    byteImgIO = io.BytesIO()
    byteImg = Image.open(('E:\\opencv\\{}.TIF'.format(a)))
    byteImgIO.seek(0)
    byteImg = byteImgIO.read()

    # img = load_img(('E:\\opencv\\{}.tif'.format(a)))
    x = img_to_array(byteImg)
    x = x.reshape((1, ) + x.shape)
    i = 1
    for batch in datagen.flow(x, batch_size=16,
        save_to_dir='E:\\opencv\\',
        save_prefix= a,
        save_format='tif'):
        i += 1
        if i > 5:
            break
    a += 1
    if a > 3:
        break


texts:
I am currently working with a lot of points of which every point has to be labeled.
I am already using the "Easy Custom Labeling" tool to move the labels and create lines to connect points with labels, as soon as there is not enough space for the labels.

Now I have different problems:

It takes a lot of work to align the labels manually to reach the goal to label every feature. Is there a faster way to to it? Very often, I let one label disappear, as soon as I move another. I think in ArcGIS there was the option to convert labels to annotations and then align these annotations very close to each other. The labels in QGIS disappear as soon as they are to close together.
As soon as I go to the Composer Manager, some labels arrange differently although I choose the same scale as I used when moving the labels manually. This means that I have to arrange the labels manually in my data view, switch to the Composer Manager and refresh the Map to see, whether it finally worked. This process takes even more time.

I am happy for every hint to label more effectively!

Some maps have a much higher density of points, I just started with Italy.


texts: Looping select by Attribute height interval and exporting to individual shapefiles using ArcGIS Desktop?

texts:
I have had a look at the existing Q&As but none of them seem to be addressing this in detail.

I have a building footprint shapefile which contains a height field. I would like to group buildings according to height and store them in separate layers.
For example I would like to iterate a Select By Attribute SQL query that selects any building between 6 and 9 metres or any building <3 metres. There is thousands of building so doing it by hand is rather unpractical.

I don't mind whether ModelBuilder or Python is used to achieve the result but I am quite new to Python.


texts: Creating ArcGIS Online webapp builder Cemetery search?

texts:
I created an ArcGIS Online webapp builder app for a cemetery project. I loaded a relationship class into ArcGIS Online that relates plot polygons to a table of interred individuals. It is a 1 to many relationship. I can get the search widget to search the related table and successfully retrieve search results from the related table, HOWEVER, when I click on the results in the dropdown list that correspond to the values in the related table to have it zoom or pan to the related polygon, it tells me no results found.

Has Esri or anyone else figured out a way within the ArcGIS Online environment how to get this to function correctly? (i.e. search related table, get related results, click on result, and pan/zoom to corresponding feature).


texts: Grouping layers using PyQGIS makes them not visible on canvas?

texts:
I've tried using the information in the following posts in order to group layers in PyQGIS (Qgis 3)

Adding layer to group in layers panel using PyQGIS?
Adding layer group using PyQGIS?
https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html

I know the second post is for QGIS 2, but I assumed that by using QgsProject in place of QgsMapLayerRegistry, things might work similarly.  However all methods have the same problem.  When I put even a single layer in a group (using the python code shown below), the layer shows up in the Layers Sidebar in QGIS but it is not visible on the canvas.  The interesting thing is that when I right click on the layer in the Layer Sidebar and select Properties, the raster I want to see shows up properly in the thumbnail image at the bottom of the Style tab.  When I add the raster to QGIS the normal way (i.e. not using Layer Trees), it shows up exactly as expected.  I was wondering if anyone had any idea what might be going on.

Method 1:

    raster_layer = QgsRasterLayer(filename, name)
    QgsProject.instance().addMapLayer(raster_layer)
    root = QgsProject.instance().layerTreeRoot()
    group = root.addGroup(group_name)
    layer = root.findLayer(raster_layer.id())
    layer2 = layer.clone()
    group.addChildNode(layer2)
    root.removeChildNode(layer)

Method 2:

    raster_layer = QgsRasterLayer(filename, name)
    group = root.addGroup(group_name)
    QgsProject.instance().addMapLayer(raster_layer, False)

    layer = QgsLayerTreeLayer(raster_layer)
    group.insertChildNode(0, layer)

The only method that works to successfully add the raster to the QGIS canvas is simply:

    raster_layer = QgsRasterLayer(filename, name)
    QgsProject.instance().addMapLayer(raster_layer)

In fact even when I don't try to add the raster to a group, but try to add it to the root itself using any other method besides addMapLayer(), the raster doesn't show on the canvas.

For example:

    raster_layer = QgsRasterLayer(filename, name)
    root = QgsProject.instance().layerTreeRoot()
    QgsProject.instance().addMapLayer(raster_layer, False)
    root.addLayer(raster_layer)

Or using the example from the third link:

    raster_layer = QgsRasterLayer(filename, name)
    QgsProject.instance().addMapLayer(raster_layer, False)
    layerTree = iface.layerTreeCanvasBridge().rootGroup()
    layerTree.insertChildNode(-1, QgsLayerTreeLayer(raster_layer))


texts: Tasselled Cap Transformation on MODIS in GRASS

texts:
I would like to know how to perform a Tasselled Cap Transformation on MODIS, in GRASS.

I've found this

https://searchcode.com/codesearch/view/19063909/

but I don't know if this has been tested.

I'm a GRASS user, but if someone knows how to do it in ENVI or ArcGIS I'd be very grateful.


texts: I am preparing layout in ArcGIS 10. I have duplicated data frame to have 2 images of the same area on single page. When I enter the legend, labels are in too small digits (24). I want it to bigger than this. Even I set the digits to 24 to 36 (make it bigger) in properties. It is still the same, unreadable. So how can I make it bigger? what is the problem with it?


texts: Merged two polygons in QGIS, but can't get rid of old boundaries that remain inside it

texts:
I've merged two polygons in QGIS, but I've been unable to rid myself of the excess boundaries from the old polygons that continue to exist inside the new one. The only boundaries I would like to keep are the outermost.

The polygon in question is the green one. I've also tried using dissolve, but it hasn't helped.


texts: Exporting ActiveView to Jpeg using ArcGIS 10.2.2 fails

texts:
I have some legacy VBA code that exports the activeview to a jpeg, which needs to be updated, so it can be used in ArcGIS Desktop 10.2.2. The application uses the following code

    Option Explicit
Public Sub ExportToJpeg()
'This will export the current map as a jpeg
    Dim pMxDoc As IMxDocument
    Set pMxDoc = ThisDocument
    Dim pActiveView As IActiveView
    Set pActiveView = pMxDoc.ActiveView
    Dim strImgFolder As String

    strImgFolder = "C:\Temp"
    'If Dir(strImgFolder, vbDirectory) = "" Then
     '   MkDir strImgFolder
    'End If
    'Create an ExportJPEG object and QI the pExport interface pointer onto it.
    ' To export to a format other than JPEG, simply create a different CoClass here
    ' and change the file extension of the ExportFileName.
    Dim pExport As IExport
    Set pExport = New ExportJPEG
    pExport.ExportFileName = strImgFolder & "\" & "aaaa" & ".jpg"

    'Set the export object resolution to 96 dpi.  The default screen resolution for
    ' MS Windows is usually 96dpi, and this value is also the default resolution for all
    ' of the ArcGIS image format exporters.  The vector format exporters (PDF, AI, etc.)
    ' default to 300 dpi, so we should explicitly set resolution here to be certain we
    ' are working with a resolution of 96.
    pExport.Resolution = 96

    'For both PageLayout and Map objects, the ExportFrame property always holds a tagRECT
    ' structure appropriate for exporting.  The values in the tagRECT correspond to the width
    ' and height to export, measured in pixels with an origin in the top left corner.
    Dim exportRECT As tagRECT
    exportRECT = pActiveView.ExportFrame

    'Create a new envelope object and populate it with the values from exportRECT.
    ' We need to do this because the exporter object requires an envelope object
    ' instead of a tagRECT structure.
    Dim pPixelBoundsEnv As IEnvelope
    Set pPixelBoundsEnv = New Envelope
    pPixelBoundsEnv.PutCoords exportRECT.Left, exportRECT.Top, exportRECT.Right, exportRECT.bottom

    'Assign the envelope object to the exporter object's PixelBounds property.  The exporter object
    ' will use these dimensions when allocating memory for the export file.
    pExport.PixelBounds = pPixelBoundsEnv

    'Initialize the exporter object and store it's device context in the hDC variable.  At this method
    ' call, the exporter object will create an empty file and allocate memory based on resolution,
    ' bit depth, and pixel bounds.
    Dim hDC    As Long
    hDC = pExport.StartExporting

    'Redraw the active view, rendering it to the exporter object device context instead of the app display.
    'We pass the following values:
    ' * hDC is the device context of the exporter object.
    ' * pExport.Resolution is 96, the default resolution value for all image export formats.  Default screen
    ' resolution for MS Windows is usually 96dpi.
    ' * exportRECT is the tagRECT structure that describes the dimensions of the view that will be rendered.
    ' The values in exportRECT should match those held in the exporter object's PixelBounds property.
    pActiveView.Output hDC, pExport.Resolution, exportRECT, Nothing, Nothing

    'Finish writing the export file and cleanup any intermediate files.
    pExport.FinishExporting
    pExport.Cleanup

End Sub

Which I got from here http://edndoc.esri.com/arcobjects/9.0/Samples/Cartography/Printing_and_Exporting/ExportActiveView.htm , it works perfectly in Desktop 10.0 but not in 10.2.2

Does anyone know how to get the above code working in 10.2.2?

EDIT -
Should have mentioned, no error occurs, The code runs to the end, a JPEG is exported but it is blank with no data, even though there is data within the map.


texts: GeoTools: java.util.{List, Map} in SimpleFeature

texts:
I'd like to add a List(Map) as Attribute to a SimpleFeature when using the PostGIS datastore:

SimpleFeatureTypeBuilder simpleFeatureTypeBuilder = new SimpleFeatureTypeBuilder();
simpleFeatureTypeBuilder.add("list", List.class);

this causes:

Dec 05, 2019 8:48:51 AM org.geotools.jdbc.JDBCFeatureSource buildFeatureType
WARNING: Could not find mapping for 'list', ignoring the column and setting the feature type read only

Therefore my question is how to I register the mapping properly?


texts: How do you make map styles NOT appear in a legend?

texts:
I am using GeoServer to serve up OSM data. I am making use of the GetLegendGraphic command, but I do not want some of my styles to appear in the legend (for example, the styles that are ONLY for text labels should not get a line).

Is there any way to specify that a style NOT appear in a legend?


texts: Display GeoPandas on rasterio doesn't work inside loop

texts:
I have written a loop that suppose to open raster from file, display it as rgb and then plot on top of it GeoPandas plot.
The problem is that the plots do not appear on to of the rasters inside the loop:

shapes = gpd.read_file(r'shape/shapes.shp')
shapes.dropna(subset=['geometry'],inplace=True)

...
#opening the files with the tif
files = glob.glob(str(folder_) + '/**/*.tiff', recursive=True)

for i in files:
    the_file=str(i)
    match=re.search(r'\d{4}-\d{2}-\d{2}',the_file)
    d=match.group(0)

    #open image

    src = rasterio.open(i)
    print(src.crs)
    #create rgb visualization
    fig, ax = plt.subplots(figsize=(20,10))
    ax=show(src.read([4,3,2]),transform=src.transform,title=d,vmin=0,vmax=30000)
    shapes.plot(ax=ax,facecolor='none', edgecolor='red')#,color='red',legend=True)

the results:

no polygons on top of the image.
When I print outside the loop with ax.get_figure() it does plot the polygons:

but when I add it to the script I don't get the plots when I show the image:
    src = rasterio.open(i)
    print(src.crs)
    #create rgb visualization
    cmap = ListedColormap(['red'], name='allred')
    fig, ax = plt.subplots(figsize=(20,10))
    ax=show(src.read([4,3,2]),transform=src.transform,title=d,vmin=0,vmax=30000)
    shapes.plot(ax=ax,facecolor='none', edgecolor='red')#,color='red',legend=True)
    ax.get_figure()

(result: same image with no polygons).
Why doesn't it show the polygons inside the loop? How can I solve it?
My end goal is to display the plots on top of the image inside my loop.


texts: Adding and calculating X and Y attributes in integer raster using ArcGIS Spatial Analyst?

texts:
I have an integer raster created from an input floating point (LiDAR) raster by multiplying it by the (elevation) value of 100 (any other value higher would not work)

Thus, I can add a "double" attribute called elevation to this raster and calculate my real-world elevation by dividing the integer value field by 100. This gives me elevation attribute with 2 decimal places in this integer raster.

I would now like to add 2 more fields to this raster, namely X and Y. Unfortunately, calculate geometry is not an option.  I can't think of any other way other than converting the raster to points then converting back to raster which can take a very long time with full resolution LiDAR rasters.

Ultimately I want to end up with an integer raster that also has X,Y,Z ("double") attributes for each cell.

Is this possible?


texts: Data type for point coordinates input in PyQGIS

texts:
I am working in QGIS console. I'd like to create new point in existing layer and using two floats as coordinates [1604288.165, 642862.0068].
feat.setGeometry(QgsGeometry.fromPointXY(coords[0], coords[1]))

But getting TypeError ((
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 80, in <module>
TypeError: QgsGeometry.fromPointXY(): argument 1 has unexpected type 'float'

Into which type should I convert my coordinates?


texts: How to add a webpage link to google map marker description

texts:
I have multiple pin points on my google maps, with a small description.
I want to add a link page (on my computer) to the description and make it linkable. Whenever I clicked, a new tab opens? Is it possible?
I tried this on google map marker

description

but it just displays it as a text.


texts: Arcgis9.3 Problem with search cursor and update cursor in two different table

texts:
I am making a python code for ArcGIS 9.3 :

-I have one folder with a lot of shapefile polygon inside
   - Ech polygon shapefile have a field "Average" and one ID
   - I have a table with each field corresponding to the ID of shapefile

I want to open each shapefile, take the "Average" and put it in the good field. So for shapefil 1, I take all the averages and i put them in line 1 of the table.

To do it I ha one SearchCursor in my shapefile and one UpdateCursor in my table. But it doesn't work and it seems to come from the update function.

The code is as following and the error come after "Marqueur 6.1" which is link to the SetValue method and to the update cursor.

The mesage error is :
 Marqueur 6.1
: ERROR 999999: Error executing function.
Failed to execute (OpeTable2).

for infile in glob.glob(os.path.join(path, '*.shp')):

    MyShape = os.path.join(path, infile)
    Name = os.path.basename(MyShape)
    fileName, fileExtension = os.path.splitext(Name)

    gp.AddMessage(Name)

    if TestFieldCreate == False :
        TestFieldCreate = True
        rowsP = gp.UpdateCursor(MyShape)
        rowP = rowsP.Next()
        #row est une ligne et non une colonne
        while rowP:
            ID_BVs = rowP.GetValue(Field_ID_BV)
            ID_BV_STR = str(ID_BVs)   #Converti en string
            NameBV = "BV"+ ID_BV_STR
            rowsP.updateRow(rowP)
            rowP = rowsP.Next()
            gp.addfield (MyTable, NameBV, "double")

        #On place le curseur dans le shp

    rowsInShape = gp.SearchCursor(MyShape)
    rowInShape = rowsInShape.Next()

    while rowInShape:

        countStr =str(count)
        gp.AddMessage( "Iteration N°" + countStr)
        Moyenne = rowInShape.GetValue(Field_Moyenne)

        if Moyenne != 0 :
            gp.AddMessage("Marqueur 6")

            # On récupère les valeur pour AfoisP et le num de BV que l'on converti de la même manière que le nom des champs
            BV = rowInShape.GetValue(Field_ID_BV)
            BV_STR = str(BV)   #Converti en string
            FieldName = "BV"+ BV_STR
            gp.AddMessage("Marqueur 6.1")
            rowInTable.SetValue(FieldName, Moyenne)
            rowInTable.SetValue("NomFich", fileName)
            rowsInTable.updateRow(rowInTable)
            gp.AddMessage("Marqueur 7")

        gp.AddMessage("Marqueur 7.1")
        rowInSHape = rowsInShape.Next()
        count=count+1

    gp.AddMessage("Marqueur 7.2")
    rowsInTable.updateRow(rowInTable)
    rowInTable = rowsInTable.Next()
    gp.AddMessage("Marqueur 8")
    NumStr = str(Num)
    gp.AddMessage(Name + " Fichier N° " + NumStr)
    Num = Num +1
    gp.AddMessage("Marqueur 9")

Thanks a lot


texts: Anyone recognize this (reference-system related?) pattern of distortion?

texts:
I'm new to GIS and was making some progress until getting rather stumped on this one.  I've imported what I believe to be lat/lon data based on US FIPS county codes in weather alerts (NOAA SPC).  When I use the geodjango "simplify" method on my geometries, they render on a map according to my expectations.  But when I attempt to perform relationship queries on the raw geometries, I get distance and overlap values that do not correspond to the real world geography.  For example, the two upper blobs in this top image should be adjacent and much further south, and the lower blob should fit Florida better (the distances and overlaps in this image match my ST_Distance and ST_Overlaps values):

Here is a rendering of the "simplify" geometry, which shows how these spatial entities are truly related in the real world (pink regions):

Does this distortion look familiar to anyone?  Thanks!!


texts: Issue with mapping MODIS SIN grid to EASE 2.0 with gdalwarp

texts:
I am trying to merge and reproject MODIS 250m NDVI product (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/myd13a3_v006) to EASE2.0 grid. I am interested in projecting these to 9km EASE2.0 grid (defined by NSIDC for SMAP mission), but I have tried with the 36km EASE2.0 grid as well and I have the same problem as below.

I first use pymodis library to merge all the tiles from each month of MODIS data and generate a mosaic version of the data.
Then I use gdalwarp with the following command to project to EASE2.0 grid:

gdalwarp -overwrite -tr 9008.055210146 9008.055210146 -t_srs '+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m' -r average mosaic_input.tif output.tif

The issue is that the output.tif file does not have the same lat/lon in each pixel as the EASE2.0 9km grid. The number of columns in the longitudinal direction are the same, but the corresponding longitudes are not the same.
These are the longitudes of the center of pixels from columns 1 to 10:

179.895906749369    179.989345967666    -179.917214814038   -179.823775595741   -179.730336377444   -179.636897159148   -179.543457940851   -179.450018722554   -179.356579504258   -179.263140285961

The same issue happens with latitudes. The upper latitude in the output.tif file is 86.9278324335268 (center of the pixel), while the 9km EASE2.0 grid has an upper latitude of 84.6564407348633.

I should add that I am sure the merging with pymodis is working fine because when I use the mosaic file with gdalwarp and reproject it to EPSG 4326 (rectangular lat/lon grid) the output file has the the correct lat lon info.


texts: Using Global Variables in Python Toolbox?

texts:
How do I create a global variable in a python toolbox?

The eventual purpose is that I would like to store information related to the execution of my script across functions, but not as "parameters".  They might be boolean, or dictionaries, but in my example below, I am using a "String for simplicity.

Add text to the Something parameter, move to Something Else, and click Ok.
On execution, the value of self.whatever is added to the messages in ArcCatalog.  It's the initial value of "ABC".

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]

class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Tool"
        self.description = ""
        self.canRunInBackground = False
        self.whatever = "ABC"

    def getParameterInfo(self):
        """Define parameter definitions"""

        param0 = arcpy.Parameter(
            displayName="Something",
            name="something",
            datatype="GPString",
            parameterType="Required",
            direction="Input"
        )

        params = [param0]

        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""

        self.whatever = "XYZ"

        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""

        arcpy.AddMessage("The value whatever is: " + self.whatever)

        return


texts: ogr2ogr when importing from file geodatabase and using in Python?

texts:
I want import data from my GDB to PostGIS. It works fine on some datasets (all records -data- is imported or nothing is imported) from a gdb but just doesn't import others. I can't find any documentation on why from sites like http://www.bostongis.com/PrinterFriendly.aspx?content_name=ogr_cheatsheet

So I have
ogr2ogr -f "PostgreSQL" -a_srs "EPSG:4283" PG:"dbname=tp_test user=postgres port=54321 password=XXXXXXXXXXX host=localhost" d:/scratch/arcgis/scratch.gdb

It only seems to report errors but if it works then you just get the command line again. When I check the postgis database, some of the datasets from the gdb are imported and others aren't. They are all in the same EPSG.

Any ideas?

This is part of a larger script/tool that I have already written to read my mxd's, export and georeference the data to a single format, save the style layer, convert it to SLD and ideally the next step is to save it to a postGIS database. This can then be used in Geoserver with the correct ESRI symbology. See Converting ArcGIS style (symbology) to font for MapInfo/OpenGeo?


texts: Setting both spatial extent and resolution when using gdalwarp?

texts:
I would like to overlay two tiff image that orginally have different projections. One is MODIS sinusoidal and the other one is WGS84. I use gdalwarp to re-project the latter into the same projection as the former. After that the two images matches in size but there is still some shift.
The gdalwarp command I used is :

gdalwarp -t_srs modis.prj -tr 463.312716 -463.312716 in.tif out.tif
modis.prj contains the WKT that configures the target projection.

I think the mismtach is becasue the spatial extent (i.e., origin) between the two data does not match perfectly. I know using -te xmin ymin xmax ymax in gdalwarp can set the spatial extent. My confusion is that I think given the resolution, only one pair of (xmax,ymax) or (xmin,ymin) should be needed, no ?

My hope is after setting exactly the (xmax,ymax) in gdalwarp as the former data set, with exactly the same resolution and the same projection system, then these two data sets should perfectly match.

Does anyone have any comments/suggestions?


texts: OpenLayers 6: GeoJSON polygons not displaying without any error

texts:
I have two GeoJSON files, the data should be Cartesian already. One GeoJSON (https://jsonblob.com/a4e3c8c9-0970-11ea-98b2-a57dd05dc415) file works fine but the polygons from the second one (https://jsonblob.com/d0ada7b7-0970-11ea-98b2-afb1520d8a49) don't show up and I don't get why.

OpenLayers do not throw any errors at all so it's difficult for me to find out the issue. Btw. I don't need any map tiles etc. I just want to see my polygon map. I can see all polygons in application like QGIS or mapsharper.org.

Here is a little bit code:

var vectorSource = new ol.source.Vector({
    features: (new ol.format.GeoJSON).readFeatures(geoJsonFeatures)
});

var vectorLayer = new ol.layer.Vector({
    source: vectorSource,
    style: stylePoligon
});

map = new ol.Map({
    //renderer: (['webgl', 'canvas']),
    target: 'map',
    layers: [
        vectorLayer
    ],
    view: new ol.View({
        //zoom: 14
    }),
    interactions: interactions,
    controls: ol.control.defaults({
        attribution: false,
        zoom: false,
    })
});

map.getView().fit(vectorSource.getExtent(), map.getSize());


texts: SQL Server geography.Intersects Yields Unexpected Results (arch shape)

texts:
I've coupled ESRI's ArcGIS with SQL Server 2008 R2 in an attempt to allow users to search for geographical points using a polygon (with ArcGIS' drawing toolbar).  I've run up against an issue that's not making any sense to me.  It seems easiest to illustrate it visually:

I draw a rectangle along the top of the northern hemisphere and get results I wouldn't expect (I know there are more points up there):
I draw a smaller square within that rectangle and get the results I was expecting:

Can anyone explain why this might be happening and if there's some way around it?  It seems like a SRID mismatch, but I'm using 4326 for ArcGIS and for SQL Server.

If I switch to using geometry instead of geography the rectangle works fine.  I assume this is because geometry isn't taking into account the curvature of the earth, which isn't present on the flattened-out map.  geometry has problems of its own though (my map wraps infinitely left/right, and so some of my polygons refer to points outside 180 degrees, which geometry doesn't seem to support).

For reference, here are the results in SQL Server:

For the rectangle:

For the square:

And UnitsSearchByLocation looks like this:

ALTER PROCEDURE [dbo].[proc_UnitsSearchByLocation]
    @polyText varchar(4096),
    @wkid int
AS
    declare @g geometry
    set @g = @polyText

    declare @poly geography

    begin try
        -- try the polygon as provided, although it has the potential
        -- to be oriented incorrectly
        set @poly = geography::STPolyFromText(@polyText, 4326)
        select * from Unit where Location.STIntersects(@poly) = 1
    end try
    begin catch
        -- assume the error was a result of the orientation being wrong
        -- and attempt to correct it
        declare @validPolyText varchar(4096)
        set @validPolyText = @g.MakeValid().STUnion(@g.STStartPoint()).ToString()
        set @poly = geography::STPolyFromText(@validPolyText, @wkid)
        select * from Unit where Location.STIntersects(@poly) = 1
    end catch
GO

Although it seems irrelevant to me, I'm using this Imagery from ESRI (:
 - https://server.arcgisonline.com/ArcGIS/rest/services/ESRI_Imagery_World_2D/MapServer

Thanks for any help you can provide..


texts: Only one messagebar created but says there's "1 more"

texts:
Using QGIS 3.14.1. If I type
iface.messageBar().pushWidget(iface.messageBar().createMessage( u'Hello World' ), Qgis.Info)

into the python console, it says there's "1 more message". But no other messages are created.
Is this a bug?


texts: Creating End Nodes using ArcGIS for Desktop?

texts: Is there any way, other than Planarizing the lines, to create an end node automatically when snapping one line to the vertex of edge of another line.


texts: How to read the bounding box of an osm file without loading it?

texts:
I noticed a pattern in .osm files : the third line tends to be its bounding box. Typically, the xml element that contains the bbox is called bounds.

I'd like to read this element without loading the whole file. Indeed, loading 80GB of data only to read one row seems pretty inefficient.

Is there a tool that lets you do this simple query ?


texts:
I'm using cartodb.js in my HTML file. I have this code:

var layerSource = {
            user_name: 'myuser',
            type: 'cartodb',
             sublayers: [{
             sql: "SELECT * FROM post",
               cartocss: '#post {marker-fill: red;}'
                        }]
                   }

Does someone know how I could get each single markers to be draggable and editable?

Now I am sure the Google Maps API can display my base map but with this library, I can't select the markers in the CartoDB Layer.


texts:
I have a rainfall time series for 15 hours in increments of 0.250s, I need to convert it from a .txt-file to a raster.
I've tried using r.mapcalc in GRASS, however, have not been successful. I've also tried creating a temporal data set in GRASS, but the layer creation always ends up empty.
The output format needs to be a raster either asc or tif format, essentially a hyetograph in raster form.
I don't mind if this is in either QGIS or GRASS. I've tried adding the layer extent to the txt file i.e. ncols, nrows and so on, cause the rainfall is uniform over a given location but I still haven't been successful.
The format is as follows:
time,rain
0,0
0.25,0
0.5,0.06100329
0.75,0.067136294
1,0.073993471
1.25,0.081679436
1.5,0.090317861

Are there any other methods I could try out for the conversion?


texts: Determining Flow Direction in ArcGIS for Desktop?

texts:
I have a problem with the Flow Direction tool ArcGIS 10.2.2.
After running the tool on a digital elevation model in .tif format 16 bit unsigned integer.
However the result is an image of 8 bit unsigned integer without flow values. See figure.

Can someone help me?


texts:
I have a simple editing application with the following code:

startEditing : function(template) {

    var drawingTool = template.template.drawingTool;

    switch(drawingTool) {
        case FeatureTemplate.TOOL_POINT:
            drawingTool = Draw.POINT;
            break;
        }

    on(this.drawingToolbar, "draw-end", lang.hitch(this, this.createFeature, template));

    this.drawingToolbar.activate(drawingTool);
}

and

createFeature : function(template, evt) {

    var featureLayer = template.featureLayer;
    template = template.template;

    var prototype = template.prototype;

    var geometry = evt.geometry;

    var graphic = new Graphic(prototype.toJson());
    graphic.setGeometry(geometry);

    this.initAttributes(graphic, featureLayer).then(function() {

        var features = [graphic];

        featureLayer.applyEdits(features).then(function(addResults) {
            var objectIds = array.map(addResults, function(addResult) {
                return addResult.objectId;
            });

            var q = new Query();
            q.objectIds = objectIds;

            featureLayer.selectFeatures(q).then(function(features) {
                main.openForm(features);
            });
        });

    });

}

Using a point tool with the drawing toolbar, I am finding that when I click a point on the map, the createFeature function is called twice resulting in 2 features being created with the same geometry. How can I stop the draw toolbar from creating 2 features?


texts: Why is arcpy.mapping.ExportToPDF losing layer symbology?

texts:
When exporting a map it does not export the symbology within the data frame.  I am using the code snippet below.

>>> myMap = arcpy.mapping.MapDocument("C:/USERS/NISUSER/DESKTOP/mxd_map.mxd")
>>> arcpy.mapping.ExportToPDF(myMap, "C:/AgTech/agtech.pdf")

The first image is what is located in the data frame:

I am expecting the 2nd image to be identical to the 1st, but as a PDF.  The symbology changes in the process somehow.


texts: Work around for “Couldn’t load PyQt” error with QGIS 3.0.3 on Ubuntu 18.04?

texts:
I am working with a fresh install of Ubuntu 18.04 and installed QGIS using https://qgis.org/debian. I get the following error when I launch QGIS:

Couldn't load PyQt.
Python support will be disabled.
Traceback (most recent call last):
  File "", line 1, in
ImportError: /usr/lib/x86_64-linux-gnu/libQt5Core.so.5: version `Qt_5.10' not found (required by /usr/local/lib/python3.6/dist-packages/PyQt5/QtCore.so)
Python version:
3.6.5 (default, Apr  1 2018, 05:46:30)
[GCC 7.3.0]

QGIS version:
3.0.3-Girona 'Girona', 8a899c8

Python path:
['/usr/share/qgis/python', '/home/ned/.local/share/QGIS/QGIS3/profiles/default/python', '/home/ned/.local/share/QGIS/QGIS3/profiles/default/python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python36.zip', '/usr/lib/python3.6', '/usr/lib/python3.6/lib-dynload', '/usr/local/lib/python3.6/dist-packages', '/usr/lib/python3/dist-packages']

A similar error has been reported as a bug (https://issues.qgis.org/issues/19040) but I haven’t seen a workaround.

I tried to install 2.18 (ltr) but apt-get failed with:

qgis : Depends: python-qgis (= 1:2.18.20+28bionic) but 1:3.1.0+git20180529+ad1fef9+28bionic is to be installed

Any suggestions on how to get QGIS working on Ubuntu 18.04?


texts: GeoJSON add and format line features to a Leaflet map?

texts:
I'm putting some line features on a Leaflet web map.

Thus far I have been able to add point features, following this tutorial, which has me turning the GeoJSON into a JavaScript variable then kicking out the following JavaScript code:

    var skull = L.icon({
    iconUrl: 'images/blackSkull.svg',
    iconSize: [38, 95],
    popupAnchor: [1, -20],});

function fatalitiesLayer (feature,
layer) {
    layer.bindPopup("
Cyclist Fatality

name: "+feature.properties.name+
    "
date: "+feature.properties.date+
    "
location: "+feature.properties.location+
    "
cause: "+feature.properties.cause+
    "
"
    );
    layer.setIcon(skull);
}

L.geoJson(fatalities,{
    onEachFeature: fatalitiesLayer
}).addTo(map);

var linesFeatureLayer = L.geoJson(segments);
linesFeatureLayer.addTo(map);

function styleLines(segments) {
    return {
                color: black,
                weight: 10,
                opacity: .7,
                dashArray: '20,15',
                lineJoin: 'round'
            };
}

L.geoJSON(segments,{
style: styleLines
}).addTo(map);

Note that the fatalities are points and the segments are lines. As things stand the points show up on the map but not the lines.

Additionally this is what my index file looks like:

  <meta charset="utf-8">
  <title>LeafletMap</title>
  <link href="style.css" type="text/css" rel=stylesheet>

Thus far all the tutorials I've found that made much sense to me were about creating one or two line features, but I've got about two hundred in a GeoJSON file and I'm not sure how to proceed.


texts: Best way for rendering geotools in web service

texts: I'm using geotools in my web service project. I always thought geotools is somehow compatible with openlayers, but now when it actually came to implementing this functionality, i can't find even one example. Any other ideas, tools to do it? I've read that geoserver is using openlayers with geotools, but as i said i can't find any example,documentation.


texts: Creating the node ids：Is this SQL error related to geometry type？

texts:
I loaded a multilinestring layer into PostGIS Manager：

and use the SQL query, trying to get the table with the node ids：

ALTER TABLE split_4 ADD COLUMN source integer;
ALTER TABLE split_4 ADD COLUMN target integer;
SELECT assign_vertex_id('split_4', 800, 'geom', 'gid');

and got this error：

Error: Function StartPoint (geometry) does not exist
LINE 1: SELECT gid AS id, StartPoint (geom) AS source, EndPoint (geom) ...
                      ^
HINT: did not comply with the specified name and parameter types of function. You may need to add explicit type conversion.
QUERY: SELECT gid AS id, StartPoint (geom) AS source, EndPoint (geom) as target FROM split_4
CONTEXT: PL / pgSQL function "assign_vertex_id" line 23 at target EXECUTE statement FOR

Is this error related to the multilinestring geometry？

I work under WIN7-64bits, with postgres 8.4, pgadmin 3, PostGIS 2.0 installed, and download Windows Binaries of pgrouting：pgRouting-1.03_pg-8.4.2.zip.


texts: Merge multiple polygon shapefile to a single polygon shapefile

texts: I'm looking to merge all internal polygons in a gridded shapefile to one solid polygon shapefile while retaining the external shape.


texts: I am trying to join a CSV table to a shapefile in QGIS 2.12.0. I have exported a table from Excel to CSV format, but when opening the table I noticed that my unique ID field in Excel is in text format, is changed to number in QGIS. Since my shapefile has that same ID field as a text string I cannot join them. I have tried to change the format to string but it does not let me. Does anyone know how to do this?


texts: I get a constant for the whole raster by subtracting two rasters in R and QGIS

texts:
I have two raster in brick format, as a result of previous operations, they have the same characteristics. When I try to make a simple subtraction between the two, it gives me a constant in all my study area. when I do other operations like multiply, add or divide it works perfectly, but not with subtraction. Anyone have any idea what might be happening?
reading the previous queries I have tried replacing NA values, with different raster formats and nothing. the input values of the subtraction are obtained by different functions.
Rtoa_B1 <- Rtoa_F(d, RadianceB1, Eo,Us)
Rsuperf_B1 <- Rsuperf_F(RadianceB1,Lr1,Tv,Eo,d,Oz,Tz,Edown)
Ratm <- Rtoa_B1 - Rsuperf_B1

the classes of the obtained files are
> Rtoa_B1
class      : RasterBrick
dimensions : 118, 111, 13098, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.004500776, 0.004518253  (x, y)
extent     : -75.7172, -75.21762, 5.978564, 6.511717  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +a=6378137 +rf=298.25722293287 +no_defs
source     : memory
names      : B1_MOD02HKM.A2019001.1605.061.2019002013822
min values :                                  0.08510739
max values :                                   0.3040431

> Rsuperf_B1
class      : RasterBrick
dimensions : 118, 111, 13098, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.004500776, 0.004518253  (x, y)
extent     : -75.7172, -75.21762, 5.978564, 6.511717  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +a=6378137 +rf=298.25722293287 +no_defs
source     : memory
names      : B1_MOD02HKM.A2019001.1605.061.2019002013822
min values :                                   0.0617993
max values :                                    0.280735

> Ratm
class      : RasterLayer
dimensions : 118, 111, 13098  (nrow, ncol, ncell)
resolution : 0.004500776, 0.004518253  (x, y)
extent     : -75.7172, -75.21762, 5.978564, 6.511717  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +a=6378137 +rf=298.25722293287 +no_defs
source     : memory
names      : layer
values     : 0.02330808, 0.02330808  (min, max)

the functions
Rsuperf_F <- function(radiance,Lr,Tv,Eo,d,Oz,Tz,Edown){

  return((pi*(radiance-Lr)*(d**2))/(Tv*((Eo*cos(Oz)*Tz)+Edown)))
}
Rtoa_F <- function(d,radiance,Eo,Us){

  return(  (pi*(d**2)*radiance)/(Eo*Us)   )
}

function values are previously calculated


texts: How to delete a column from the attribute table of a SpatiaLite layer?

texts:
I want to delete a column from Attribute table but the function "Delete columns" isn't highlighted. (When I am adding new layers I use New SpatiaLite layer, not Shapefile layer).

So I have installed the plugin Table Manager but when I delete a column using this Table Manager the function "Save" isn't highlighted. Only "Save as" is possible. Is there any advise how I could delete a column?


texts: Creating convex hull of set of points using QGIS

texts:
I have a shapefile layer with a set of points.

I would like to enclose them (or turn them) into the polygon (with the line in the middle).

I tried this:

Converting cluster of Points to Polygons using QGIS?

Creating Minimum Convex Polygon - Home Range from Points in QGIS

which led me to the Convex Hull tool.

Unfortunately, it looks like, that this operation can be applied to the cluster of points instead of single located points.