texts: How to add on click information on places in any layers in OpenGeo Suit 4.5?

texts: I recently installed OpenGeoSuite 4.5 and I am very new in to it. I uploaded my .shp file in it and it showing that layer currently. Now I want to add some information(name and a little info etc) of those places . Which i want to show by clicking on that particular place like state or  building etc. Can I also add markers and images on my map? So please give me steps to do it.


texts: How to install postgis extension using an existent version on postgresql on OSX 10.9.4

texts:
I'm trying to install postgis using my existent version of postgresql.

I tried doing: brew install postgis. But I'm getting the following error trying run create extension postgis; query

PG::UndefinedFile: ERROR:  could not open extension control file "/opt/boxen/homebrew/Cellar/postgresql/9.3.2-boxen/share/postgresql/extension/postgis.control": No such file or directory

This is happening because brew is installing a fresh version of postgresql. When I inspect on the new postgres path I can find the missing postgis.control file.

So, the question is: how can I force postgis to use my existent postgresql installation?

(I'm using boxen, so It's not an option remove my actual postgres installation)


texts: Analysing landuse data in a buffer layer

texts: In my project I want to analyse the landuse around some points within a certain radius. So I created a buffer around my points but I have no idea how QGIS could probably tell me what different landuses are laying inside of the buffers (like 2 km² forest, 3 km² residential etc.). I have already loaded a vector layer which shows the different landuses in the whole area in general on the map.


texts:
I'm working on a custom toolbar for an extension which will have operations that work upon a specific, user-selectable layer.

Is it possible to directly use the layer selector used by ESRI's toolbars in your own toolbar?

For example, the 3D Effects toolbar contains a "Layer:.." combo box for selecting the layer on which to operate.  It includes a full tree of the scene, filtered to certain types of layers.  Is this control a reusable object?


texts: Generating single GPX File from two feature classes in ArcGIS Pro

texts:
I currently have two feature classes. One is a point feature class consisting of sampling stations for a citizen science project. The second is a line feature class consisting of the route that an observer would take to get to those points. The two  features correspond to each other with the route number identifier (could work as a join field) but since they're two different types of feature classes I'm not sure how they would be combined.

The end goal would be to generate a geospatial information package for each route, which would consist of the points and route for each route and only for that route. It'd have to be 100 or so files in an easily accessible filetype (preferably .gpx but .kml/.kmz would work too). How could I go about this?

For reference, I'm using ArcGIS Pro and don't really have access to the licensed features.


texts: How do I transform custom coordinates to pixel in vector layer in OpenLayers 6

texts:
I have a XYZ layer for tiles(generated from an 1065 * 950 Image using gdal2tiles), and a vector layer for features, I want to use my own pixel transformation to set the coordinates of the features.
 this.extent = [0, 0, this.imgWidth, this.imgHeight];

 this.projection = new Projection({
  code: "warehouse-tile-image",
  units: "pixel",
  extent: this.extent,
});

this.view = new View({
  projection: this.projection,
  center: getCenter(
   this.extent
  ),
  zoom: 2,
  maxZoom: 3,
});
this.tileSource = new XYZ({
  url: "assets/map/{z}/{x}/{y}.png",
  minZoom: 0,
  projection: this.projection,
  maxZoom: 3,
  tileSize: 256,
});
 this.layer = new VectorLayer({
  source: new VectorSource(),
});

I want to use this d3 scale function to set the pixel coordinates:
//The coordinate of data is in meter
const scale= 0.05
const max_x = scale* this.imgWidth + this.origin.x;
const max_y = scale* this.imgHeight + this.origin.y;
this.xScale = scaleLinear()
  .domain([this.origin.x, max_x])
  .range([0, this.imgWidth]);
this.yScale = scaleLinear()
  .domain([this.origin.y, max_y])
  .range([this.imgHeight, 0]);

So I can transform the custom coordinates to pixel by using

this.xScale(x) = pixel x

and

this.yScale(y) = pixel y

This is how I do in Leaflet in order to re-projection:
    const denominator = Math.pow(2, this.maxZoom);
    L.CRS.MySimple = L.Util.extend({}, L.CRS.Simple, {
        /**
         * max zoom level >= log2(max(width,height)/tilesize)
         * a,c = 1 / 2^max zoom level
         **/
        transformation: new L.Transformation(
            1 / denominator, //a
            0, //b
            1 / denominator, //c
            0 //d
        ),
    });
 this.mapOptions = {
        layers: [
            L.tileLayer(this.tileUrl, {
                maxZoom: this.maxZoom,
                tileSize: 256,
                noWrap: true,
                bounds: bounds,
            }),
        ],
        zoom: this.maxZoom,
        crs: L.CRS.MySimple,
        maxZoom: this.maxZoom,
        minZoom: 0,
        maxBounds: bounds,
        zoomControl: false,
        center: L.latLng([this.imgWidth, this.imgHeight]),
    };
 private transformCoord(coordinate: { x: number; y: number }): number[] {
    //This function will transform custom coordinates to pixel coordinate using d3 scale
    return [this.xScale(coordinate.x), this.yScale(coordinate.y)];
}
//Feature geometry
geometry: {
                type: "Point",
                coordinates: this.transformCoord(data.pose.position),
            },

And then I just use those scale function to transform the coordinates to pixel coordinates and feed to the data layer(Leaflet-realtime which use a Geojson layer).
However, it is different in OpenLayers which seems to do the pixel transformation automatically for me.
My problem is I don't understand how the projections work in OpenLayers, I have no idea how to do the same as in Leaflet.


texts: Calculating values for new field using values from another field as values in dictionary?

texts:
I am using the ArcMap application of ArcGIS Desktop 10.5. Im trying to use field calculator to return a max value between other fields in the feature based upon a condition. I am not new to python, but I am new to implementing it in ArcMap.  Using field calculator, it seems as though you can reference that row's value for a particular field by delimiting the field names with '!'. What I would like to do is take the numeric amounts for how many people speak a given language per tract feature and then return a string of the most common language name as long as the language is spoken by at least 5% of the population. An example of a python code that I assumed would work is first having a dictionary with keys as the text string for the language name and dictionary values by referring to that field name with '!'. Then putting this dictionary and a function in the 'Pre-logic Script Code' and then in the field calculator passing in this dictionary into the function to return a string.

Pre-logic Script Code:

Totals = {'Total':!Estimate__Total_!, 'Total only English':!Estimate__Total____Speak_only_English!, 'Spanish': !Estimate__Total____Spanish_!, 'French Haitian or Cajun' : !Estimate__Total____French__Haitian__or_Cajun_!, 'German' : !Estimate__Total____German_or_other_West_Germanic_languages_!, 'Russian Polish or Other Slavic Languages' : !Estimate__Total____Russian__Polish__or_other_Slavic_languages_!,'Other Indo European Languages': !Estimate__Total____Other_Indo_European_languages_!,'Korean': !Estimate__Total____Korean_!,'Chinese': !Estimate__Total____Chinese__incl__Mandarin__Cantonese__!, 'Vietnames': !Estimate__Total____Vietnamese_!, 'Tagalog':  !Estimate__Total____Tagalog__incl__Filipino__!, 'Asian': !Estimate__Total____Other_Asian_and_Pacific_Island_languages_!, 'Arabic': !Estimate__Total____Arabic_!, 'Other': !Estimate__Total____Other_and_unspecified_languages_!}

 def mostCommon (dic):
    most = 0
    lang = ''
    tot = dic['Total']
    for k,v in dic.items():
        if (v > most) and (v >= (tot *.05)):
            most = v
            lang = k
    if lang == '':
        return 'None'
    else:
        return lang

and then in my field calculator box below I simply pass the dictionary into the function

mostCommon(Totals)

However, this gives me a parameter error.

Is there any other way for this to be done in field calculator?

I have also researched using UpdateCursor, but im not quite sure how I can reference a value in a particular field for the row and add it to a dictionary as I have above.


texts: Is there a limit to the number of polygon points you can export to CSV in QGIS 2.4?

texts:
Using QGIS 2.4

I am currently importing a csv file with ~30,000 polygons (polygons in WKT format). I also have a United States shapefile. I have used the Clip function to get only the polygons that are in the US. When I go to save the new clipped shapefile as a CSV file with GEOMETRY set to AS_WKT, sometimes the number of lat/lon points is very long and so it cuts off, not finishing the polygon or multipolygon. Is there a way to set the export to include all points?


texts: Are there any resources or standards when writing a GIS data specification?

texts:
I have been searching for sources of GIS data specifications to assist with documenting a number of datasets that I'm currently working on. I'm used to looking at Ordnance Survey specifications, such as the one located here: http://www.ordnancesurvey.co.uk/docs/technical-specifications/os-mastermap-water-layer-technical-specification.pdf.

I'm guessing that a good specification will be one which is independent of applications and formats, and describes the data in the detail which is necessary for its purpose and intended audience.

Are there any resources to assist with the creation of GIS data specifications? Is there an industry standard way to do them so that they can be easily understood by both users and GIS developers? I'd also be interested to know what kind of documentation people do generally create for datasets which are to be widely used within their organisations.

Many thanks.


texts: Unable to create shapefile using shapefile.py

texts:
I am trying to create shapefile from csv data using shapefile.py but i am constantly getting this error,

Traceback (most recent call last):
  File "createShp_Davis.py", line 31, in <module>
    w.save(File_loc)
  File "/usr/lib/python2.7/dist-packages/shapefile.py", line 862, in save
    self.saveShp(target)
  File "/usr/lib/python2.7/dist-packages/shapefile.py", line 829, in saveShp
    self.__shapefileHeader(self.shp, headerType='shp')
  File "/usr/lib/python2.7/dist-packages/shapefile.py", line 586, in __shapefileHeader
    raise ShapefileException("Failed to write shapefile bounding box. Floats required.")
shapefile.ShapefileException: Failed to write shapefile bounding box. Floats required.

For some reason it is not saving the shapefile properly, it is able to create .shp file but that's it. There are no .dbf or .shx. I believe i am making some mistake somewhere in the code, but not sure where. Here is my code.

#!/usr/bin/python
import csv
import shapefile
import BIL

File_loc = raw_input("Enter the file name and folder location: ")

w = shapefile.Writer(shapefile.POINT)
w.field('Species_Id')
w.field('Species')
w.field('MeanAnnTemp')
w.field('TempSeasonality')
w.field('MinTempColdestMonth')
w.field('AnnPrecipitation')
#w.field('Elevation')
w.field('Latitude')
w.field('Longitude')

with open('treeMap.csv', 'rb') as f:
    reader = csv.reader(f)
    for rec in reader:
        w.point(rec[6],rec[5])
        lat=[];lon=[]
        lat.append(rec[5]);lon.append(rec[6])
        lat=[float(i) for i in lat];lon=[float(i) for i in lon]
        MAT=BIL.lookupBIL(r'./worldclim/bio1',lat,lon)
        TSD=BIL.lookupBIL(r'./worldclim/bio4',lat,lon)
        AnnTMin=BIL.lookupBIL(r'./worldclim/bio6',lat,lon)
        ANNPREC=BIL.lookupBIL(r'./worldclim/bio12',lat,lon)
        w.record(rec[2], rec[3], MAT[0], TSD[0], AnnTMin[0], ANNPREC[0], rec[5], rec[6])
    w.save(File_loc)

print "Shapefile Successfully created"


texts: Converting raster to array then CSV file using ArcPy?

texts:
I have a raster file and I would like to convert it to an array and then save it as a csv file with using ArcPy.
The reason is that I want to do some statistical analysis over each row and column in the output csv file . I have around 100 rasters that I would like to have at the end 100 csv files for the further analysis. this is why I need to use arcpy in this case.
As you can see the raster output is tilted and the first row in the top part of the raster is not fully straight and it mixed with NoData as well. I would like to have the output file without the local coordinate system.
the raster resolution is 25 cm :

the example of the csv file:


texts: OpenLayers 2.13 : How to display a local tiled raster?

texts:
I am trying to display (on OpenLayers 2.13) a local tiled raster I made with gdal2tiles.py. Unfortunately it's all blank or pink. Here is my layer definition:

var tiles = new OpenLayers.Layer.XYZ("Local tiles", "tiles/${z}/${x}/${y}.png", {
    isBaseLayer: true,
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    maxExtent: [-75.1729846, 39.9439032, -75.1529085, 39.9599679]
});

map.addLayer(tiles);

What is wrong?

Here is a link to my OpenLayers test page : http://popopopopo.voila.net/index.html

And here is the original image: http://popopopopo.voila.net/raster.tif (EPSG: 3857)


texts: Integrate right to left text label expression with fields in ArcGIS 10.3

texts:
I try to insert this label expression in ArcGIS 10.3 in the layer properties. The expression has field from the attribute table and some text(right to left orientation) that i added (it isn't in some filed of the attribute table):

'אאא ' + [Muni_Heb]

I also tried to add 'u' to this expression:

u'אאא ' + [Muni_Heb]

but i get an error:

I don't understand what my mistake is.


texts: How to get the mapserver WFS OGR driver to encode unicode in geojson

texts:
Running mapserver 6.4.1, gdal 1.14.3 I try to output some WFS features using outputformat=geojson. The relevant bits of the mapfile are appended below.

Everything works fine, but unicode characters are still encoded as ISO-8859-1 when using geojson. The OGR/GML driver outputs utf-8 correctly. OGR/CSV also fails and uses ISO-8859-1.

I get for example "gatenavn": "LangÃ¸rhÃ¸gden" with outputformat=geojson and <ms:gatenavn>Langørhøgden</ms:gatenavn> in GML

Is there any dataset or layer creation option or mapfile parameter I am missing?

 (in WEB/METADATA)
 "wfs_getfeature_formatlist" "geojson,csv,ogrgml"
 "wfs_encoding" "UTF-8"

 OUTPUTFORMAT
  NAME "geojson"
  DRIVER "OGR/GEOJSON"
  MIMETYPE "application/json; subtype=geojson"
  FORMATOPTION "STORAGE=stream"
  FORMATOPTION "FORM=SIMPLE"
  FORMATOPTION "LCO:COORDINATE_PRECISION=5"
 END


texts: Extract mean for values in a raster above certain threshold in R

texts:
How to extract mean from a raster using extract function in R while considering raster cell values above 0?

I have multiple rasters with some negative cell values and no-data. I am able to ignore no data using na.rm = T. But don't know how to ignore negative cell values or apply  greater than 0 threshold.

Below is my code:

library(rgdal)
library(raster)
setwd("E:/R_data")
pointg <- readOGR("E:/R_data","euro")
rasterfiles <- list.files(".", pattern = "tif")
for(j in 1:length(rasterfiles)){
  curRaster <- raster(paste("E:/R_data", rasterfiles[j], sep ="/"))
  csvDatatoWrite <- data.frame(pointg$Name, extract(curRaster, pointg, fun=mean, na.rm=TRUE))
  names(csvDatatoWrite)[2] <- substr(rasterfiles[j],1,7)
  csvst <- csvDatatoWrite
  csvtrans <- t(csvst)
  write.table(csvtrans, file = "euro.csv", append = TRUE, sep = ",", row.names = TRUE, col.names = FALSE)
}


texts:
I am currently creating a polyline shapefile and need to split a line so as to add in a new off shoot. Would anyone happen to know how to split a polyline in QGIS?

I have tried the split features tool in the advanced digitising toolbar with the following result. It is split at the node but the rest of the line is aligned to a new bearing and the parent section I wanted keeping was lost.


texts: Finding the overall weighted centroid from other centroids

texts:
I'm new to python, and I'm trying to find the overall weighted population centroid of a state, using the unweighted centroid from counties. These unweighted centroids already have a popualtion field ("totpop") in the attribute table. The tool runs, however, it outputs the overall centroid (in pink) in the incorrect spot. It falls on the last county's centroid.

import os
import arcpy
import math
from arcpy.sa import *
from arcpy import env

with arcpy.da.SearchCursor(InPutPoints,["STATE_NAME","YCoord","XCoord","totpop", "NAME", "POLY_ID", "OBJECTID", "STATE_FIPS", "CNTY_FIPS", "FIPS", "FIPSnum","FIPS_NUMER","SHAPE@AREA","SHAPE@X", "SHAPE@Y"]) as cursor:

    Total_X1 = 0
    Total_Y1 = 0
    TotalPop1 = 0
    cntyName1 = "cnty"
    stateName1 = "statename"
    for row in cursor:
            TheStateName1 = row[0]
            thecntyName1 = row[4]
            idpoly1 = row[5]
            idobject1 = row[6]
            stateFIPS1 = row[7]
            countyFIPS1 = row[8]
            fips1 = row[9]
            fipSnum1 = row[10]
            fipsNumer1 = row[11]
            XPoint = row [13]
            YPoint = row[14]

            print TheStateName1
            print thecntyName1

            print row

            Total_X1 +=  row[2] *row[3]
            print Total_X1

            Total_Y1 +=  row[1] *row[3]
            print Total_Y1

            TotalPop1 += row[3]
            print TotalPop1

    print ""
    print "X is: " , Total_X1
    print "POP is: " , TotalPop1
    centroid_X1 = Total_X1/TotalPop1
    print "your x centroid is: ",centroid_X1
    print ""
    print thecntyName1
    print TheStateName1

########

    print ""
    print "Y is: " , Total_Y1
    print "POP is: " , TotalPop1
    centroid_Y1 = Total_Y1/TotalPop1
    print "your y centroid is: ",centroid_Y1
    print ""

    cursor01 = arcpy.da.InsertCursor(OutPutCent,["SHAPE@X", "SHAPE@Y","Xcoord","Ycoord","totpop", "NAME","STATE_NAME","POLY_ID", "OBJECTID", "STATE_FIPS", "CNTY_FIPS", "FIPS", "FIPSnum","FIPS_NUMER" ])

    cursor01.insertRow([XPoint,YPoint,centroid_X1,centroid_Y1,TotalPop1,thecntyName1,TheStateName1,idpoly1, idobject1, stateFIPS1, countyFIPS1, fips1, fipSnum1, fipsNumer1])

print" Done"


texts:
How to get a bounding box (in the form of a Polygon, or anything else where it's easy to extract the corner coordinates from) of a given Geography item? Something that would look like:

SELECT Boundary(geog) AS boundingBox FROM segments;

Where geog is a MULTILINE geography.

Seems like a trivial question, but my experience so far is limited and I did not have any luck with the docs (http://www.postgis.org/documentation/manual-1.5/reference.html ) Tried ST_Boundary and Box2D but they do not work on geographies.


texts: How to get the rectangular coordinates from list of line string

texts:
I have a linestring geometry in my PostgreSQL database. I want to get the corner coordinates of the rectangular bounding box. Example

"0102000020E610000002000000CCDA5CEB41D323406E2A9CEE86C54A40CCDA5CEB41D323406E2A9CEE86C54A40"
"0102000020E610000002000000AF5FB01BB6CD22407C7E18213C8E4A40AF5FB01BB6CD22407C7E18213C8E4A40"

Until now, I have executed the query

SELECT ST_Extent(geomtery) as bextent from warning_060219;

result

BOX(6.84285 51.37265,11.55307 54.36356)

How can I get the four corner coordinates of the rectangle?


texts: Determine If Polygon exists/intersects inside a BBOX


texts: Understanding terminologies on Esri ArcSDE

texts:
I am new to ArcSDE.

Can anyone help me understand the differences between versioning, replication, synchronization, archiving, editor tracking, and schema change?

Their concepts are similar to me.

Here are more specific questions, but you don't have to stick with answering all of these questions as I just want to understand the differences:

1) What I know is, by enabling archiving, all edits will be saved. Editor tracking allows you to see all edits. Then, why do you need both? What are the differences between archiving and editor tracking?

2) Is versioning only enable after enabling archiving or replication?

3) Is versioning like auto-save in Microsoft program? As Microsoft programs have an auto save function, you can access the file--choosing one among the various versions by saving time--when the file is suddenly closed without saving. In that case, the versioning is temporary process?

4) I get synchronization is matching a version of a file between SDEDev and SDEPub. But there are versioning and replication.. why do you need the synchronization?

5) I was told that a new replica is better to be created after schema changes. Then why do you need versioning and editor tracking?

I think I am really confused.


texts: Automatically determine polygon overlap 'islands'

texts:
I have a shapefile with many polygons. Some partly overlap with each other, some are not overlapping. I'm  considering a polygon to be overlapping given a minimum area of overlap. 'Islands' of partly overlapping polygons occur. I'm interested to find out how many islands are there and which polygons are forming them.

I can use ArcGIS 10.2 and/or OGR, both with Python but I'm not sure how to proceed. I know there is an overlap function in OGR but it just finds overlapping polygons. Is there maybe functionality already available which does the above?

Example:


texts: App for recording waypoints and photos for GIS on android

texts:
I'm looking to collect photos and notes during a field trip to add to (just for comparison) a landcover raster dataset I'm working with.

Ideally I want to visit certain field sites, take photos with my android phone, and afterwards be able to locate each photo on my map/GIS when I'm back home.

The standard location saving feature of the android camera isn't reliable enough.

I tried MapIt, which has a lot of features I don't need, but is costly. MapPad doesn't seem to permit exports into suitable formats.


texts: Creating a shapefile for footprint of raster's valid data areas with GDAL

texts:
I would like to show the borders of the areas of a given ASCII DTM file for which there is data. This problem has been discussed here previously but the recommended solution, the Image Boundary plugin for QGIS, is now considered deprecated. Furthermore, I need to be able to script it, so using the GDAL tools from the command line is probably a better plan for me.

I also have an additional restriction: I need to be able to serve this file on GeoServer, and ultimately display it using Leaflet, so a shapefile or GeoJSON could be good choices for the file format of the final output.

What command should I use?


texts: Creating surface area regions of minimum size from raster file using ArcGIS for Desktop?

texts:
I have these in a raster of Ontario province in Canada (both tiled and mosaiced versions) in ArcGIS

Also I have a list of dam sites. I created buffer zones around these sites, now I want to generate potential sites of plane areas (with low elevation distance between points of each area). Basically the output should be a shapefile with all these polygons for potential sites.

I am new to ArcGIS and GIS in general, all the examples and questions I see they already have shapefiles with sites... what emphasized text I need is to generate a shapefile with potential sites with a minimum area.


texts: Is there a way to use arbitrary, non "plugin" libraries inside Folium / JupyterLab?

texts:
It appears Folium supports certain plugins through folium.plugins
Is it possible to use an arbitrary third-party library in Folium / JupyterLab? I want to use sidebar-v2.
I know I can generate a map in Folium and tweak it outside. Looking to avoid that.


texts: How can I get an error map of an interpolation in ArcGIS?

texts:
I have about 10 point shapefiles containing a lot of points each. (I would prefer not to join them into one shapefile, it's pain). And I am using interpolation method Topo to Raster to create a DTM. This method can deal with many separate datasets. I'm using ArcGIS 9.3 (and have access to 10).

I would like to get errors of the estimates, to check if they have normal distribution.
How can I get the interpolation errors?
Like the -e flag option in GRASS.

Is it available only for kriging?
Is it the only way, to compute them manually comparing each dataset to the resulting raster? Thank you for advices.


texts: Updating certain values in a MATERIALIZED VIEW

texts:
I have created a materialized view and now want to replace a certain value in one of the columns.
How can I do that?
Perhaps I should use INSERT INTO, but I do not really know how to use this.

The road status value I want to replace is a BLANK and not a NULL value.

My code is:

UPDATE konflikt.vejstatus_dissolve
SET roadstatus = 'Not mentioned'
WHERE
   roadstatus = '' ;

Actually, the neatest thing would be if the update of data could be integrated into the materialized view.

CREATE MATERIALIZED VIEW konflikt.vejstatus_dissolve AS
       SELECT row_number() over () as gid,
       ST_UNION(the_geom)::geometry(MULTILINESTRING,25832) as the_geom,vejstistatus
       FROM_05_veje_trafik.vejstatus
       GROUP BY roadstatus;

What do I need to do?


texts: GDAL's Python wrappers don't load correct asc header informaton

texts:
I have a Python script that I am using to load .asc raster files and then work with the contents therein. However, it isn't loading the correct information that is specified in the headers. The .asc file that I am using has the following for headers:

ncols        2131
nrows        2220
xllcorner    670813.505443912582
yllcorner    5480455.007973534055
cellsize     3.000109992417

But when I run the following simple script to obtain the values of xllcorner, yllcorner, and cellsize:

from osgeo import gdal

dempath = 'My_DEM.asc'
dem_file = gdal.Open(dempath)
dem_transform = dem_file.GetGeoTransform()
print(dem_transform)

I get the following as output:

(670813.5054439126, 3.000109992417, 0.0, 5487115.2521567, 0.0, -3.000109992417)

So as we can see, the yllcorner value that is loaded from the asc file isn't 5480455.007973534055, or any rounding off of it, that I would expect it to be. Instead I am getting 5487115.2521567, which is far off enough that it causes problems for me later on. I can import the correct header information with a variety of ways that are not using the GDAL Python wrappers, but ideally I'd like to keep using the GDAL Python wrappers. Is there something special I am forgetting to specify in my script like projection information? Is this a known bug in GDAL? How can I resolve this issue without having to write my own asc file parser?


texts: PostGIS ST_RECLASS using Lookup Table not working properly

texts:
Link to question/answer similar to what I am trying to do

I am trying to run ST_RECLASS on a raster to retain certain pixel values while setting others to NoData(0).  I have a function that creates a lookup table and then I embed the lookup table within the UPDATE query using ST_RECLASS to try to reassign the output pixel valuese.

The raster, before relcass attempt, has these values: (BEC_ZONE_CLIP_RASTER_RASTER)

My lookup table contains the following values: (ZZ_TMP_LOOKUP)

When I execute the following line of code:

UPDATE BEC_ZONE_CLIP_RASTER_RASTER SET rast = ST_Reclass(rast, ROW(1, (SELECT string_agg(concat('[', val, ']:', out_val), ',')  FROM zz_tmp_lookup), '4BUI', 0)::reclassarg);

I end up with the pixels that represent 17 and 21 values reclassed, but they are set to value 15.  Can somebody help me to understand why this is happening?


texts:
i just tried the plugin Qgis2web.
The result looks already very nice

But to optimize to lay out , i have some questions.

How is it possible to display a title on the map?
How is it possible to display a logo and the name of your organisation in the right corner?

I changed already the code for the lay out in the pop up content but how is it possible to create a white space  between 'voornaam' and achternaam'.

And a last little thing: If i haven't a value in the variable 'National_1'it display 'null' instead of a white space..

var popupContent = 'Naam' + Autolinker.link(String(feature.properties['voornaam'])) +
             Autolinker.link(String(feature.properties['Achternaam'])) +
            'Geboortedatum' + Autolinker.link(String(feature.properties['Geboorteda'])) +
            'Nationaliteit' + Autolinker.link(String(feature.properties['Nationalit'])) +
             Autolinker.link(String(feature.properties['National_1'])) +
             'Positie' + Autolinker.link(String(feature.properties['Positie']))
             '';
            layer.bindPopup(popupContent);
        }


texts:
I am trying to set up a layout with the atlas-function in QGIS. I am using also the INPUT-mobile App. With the app it is possible to take pictures.  I want these pictures to be displayed using the atlas-function.
For INPUT to do make the fotos work, you have to set up the Widget type = Attachment. Also one has to use relative paths. In Integrated Document Viewer section the type is set to Image.
Now when I take a pic INPUT automatically writes in the attribute table and it looks like this:

To use the atlas function i have already discovered, that i can select the fotos from the attribute table. However, since there is not a complete path, it does not display the foto correctly. Now I am looking for a solution to fix this.

I thought the best way would be to store the atlas shape file in the same folder with the pictures and then create a path from the file directory somehow combining the information with the Foto file name.
How can I create this path or do you have a better, more convenient idea?


texts: How to upload a csv file into esri map with coordinates and plot in an esri map

texts:
I have got a sample code to plot the points on an esri map using a csv file.
But the csv file is maintained in a server, i'm unable to upload it from my local system.Please let me know how i can refer the csv file in my local system and plot the coordinates?

    Simple Map

      html, body, #map {
        height: 100%;
        width: 100%;
        margin: 0;
        padding: 0;
      }
      body {
        background-color: #FFF;
        overflow: hidden;
        font-family: "Trebuchet MS";
      }

      var map, csv;

  require([
    "esri/map",
    "esri/layers/CSVLayer",
    "esri/Color",
    "esri/symbols/SimpleMarkerSymbol",
    "esri/renderers/SimpleRenderer",
    "esri/InfoTemplate",
    "esri/urlUtils",
    "dojo/domReady!"
  ], function(
    Map, CSVLayer, Color, SimpleMarkerSymbol, SimpleRenderer, InfoTemplate, urlUtils
  ) {
   urlUtils.addProxyRule({
      proxyUrl: "/proxy",
      urlPrefix: "earthquake.usgs.gov"
    });
    map = new Map("map", {
      basemap: "gray",
      center: [ -60, -10 ],
      zoom: 4
    });
    csv = new CSVLayer("http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_week.csv", {
      copyright: "USGS.gov"
    });
    var orangeRed = new Color([238, 69, 0, 0.5]); // hex is #ff4500
    var marker = new SimpleMarkerSymbol("solid", 15, null, orangeRed);
    var renderer = new SimpleRenderer(marker);
    csv.setRenderer(renderer);
    var template = new InfoTemplate("${type}", "${place}");
    csv.setInfoTemplate(template);
    map.addLayer(csv);
  });
</script>


texts:
I am building ArcMap Add-Ins using ArcGIS 10.2 and I know that the AddIns won't work on lower versions of ArcGIS.

Is it possible to build the AddIn for ArcGIS 10.0 or 10.1 using ArcGIS 10.2?
I am using the "Python Add-In Wizard".


texts:
I am creating a simple offline navigation app in android. I have camera information(coordinates and range) and i would like to navigate the routes to avoid being caught on cameras. I have seen great example https://graphhopper.com/blog/2015/11/15/publess-routes-graphhopper-edition/ here.

What I have tried is that I have changed weights of the edges direclty connecting to cameras (problem here - is the modified weights persist even after restarting/uninstalling app).

What i need is in addition to the edges directly connecting to cameras, I would like to change weights of edges that lies in circles centring at the cameras.

I think I have read somewhere before that mapsforge provides a library with geospatial opertaions and I can't seem to find it again now.


texts: How to resize a polygon (rectangle or circle) in Qgis?

texts:
Yes I used the Rectangles Ovals plugin to create a rectangle then rotate it and move it nto position. But if you want to resize it that is where the fun starts because it is just a polygon with four corners, you can move them all independently so it is very hard to keep proportion or angles correct.

A rectangle is four lines at right angles. To make it longer or wider you just need to be able  to increase the length of two of the sides. A circle has many more "sides" so this would be harder to implement using measurements alone, but it could be done this way.

I don't know if anyone has implemented a plugin that can do this resizing. I just know dragging the nodes to new positions is very tedious especially trying to keep a circle circular. Something that let me lengthen the individual sides is a quick-and-dirty way of increasing a polygon but there may be some better option. But the Rectangles and Ovals plugin only lets me draw new ones and if the new one is the wrong size I have to start again.


texts: GetFeature WFS call to see if a particular feature exists - deliver WMS content based on the response

texts:
If I want to send a WFS query to my GeoServer server (A) from another server (B) - not any browser - what interface would I usually use? I am only familiar with accessing the GeoServer server from a browser via JavaScript and OpenLayers.

I run PHP on Server B. Browser communicates directly with server B. Server B delivers OpenLayers JavaScript to the browser. Usually, OpenLayers JavaScript communicates from the browser directly with my GeoServer server (A). Now, I want server A to send a query to server B, that obtains some data directly which would have been otherwise gone the longer way through the browser and OpenLayers.

I can do it like this:

$content = file_get_contents("http://mygeoserver?service=WFS&request=GetFeature&typeName=mytType:myName");
doSomeXMLSParseOperations($content);

but I am not sure if that is the regular solution.


texts: How to load a .osm file as a basemap tiles in openlayers 3?

texts:
I am working on an app that needs to be used offline so I have downloaded a map in an .osm file, I need to load it as a basemap. I did some research and found this code I used but it shows the map as a vector layer:

new ol.layer.Vector({
            source: new ol.source.Vector({
                format: new ol.format.OSMXML(),
                url: 'Maps/Tangier_Medina.osm',
                projection: 'EPSG:4326'
            })
})

I thought at first that specifying the format

ol.format.OSMXML()

will be enough but I need it to be a basemap with tiles.

Anyone knows how please ?


texts: Can any non-Terrago products produce geopdfs recognized by the Terrago toolbar in Acrobat?

texts: Neither ArcMap nor GlobalMapper seem to be produce geopdfs that are recognized by the Terrago toolbar. Can anything besides Terrago products produce geopdfs that the Toolbar will accept?


texts:
Got a leaflet map with some features and a button with a search function. A marker is shown at the lat,long where the feature is found. A click on the feature shows a popup.

heres the code for the feature popup:

 function pop_Alle(feature, layer) {

   var popupContent = '

BEZ	' + Autolinker.link(String(feature.properties['BEZ'])) + '

NAME	' + Autolinker.link(String(feature.properties['NAME'])) + '

TELEFON	' + Autolinker.link(String(feature.properties['TELEFON'])) + '

OEFFNUNG	' + Autolinker.link(String(feature.properties['OEFFNUNG'])) + '

WWW	' + Autolinker.link(String(feature.properties['WWW'])) + '

EMAIL	' + Autolinker.link(String(feature.properties['EMAIL'])) + '

STRASSE	' + Autolinker.link(String(feature.properties['STRASSE'])) + '

ORT	' + Autolinker.link(String(feature.properties['ORT'])) + '

BEMERKUNG	' + Autolinker.link(String(feature.properties['BEMERKUNG'])) + '

PLZ	' + Autolinker.link(String(feature.properties['PLZ'])) + '

Hkat	' + Autolinker.link(String(feature.properties['Hkat'])) + '

Nkat	' + Autolinker.link(String(feature.properties['Nkat'])) + '

';
        layer.bindPopup(popupContent);
    }

    function doStyleAlle() {
        return {
            radius: 0.0,
            fillColor: '#e31a1c',
            color: '#000000',
            weight: 0.0,
            opacity: 1.0,
            dashArray: '',
            fillOpacity: 0.0
        }
    }

    function doPointToLayerAlle(feature, latlng) {
        return L.circleMarker(latlng, doStyleAlle())
    }
    var json_AlleJSON = new L.geoJson(json_Alle, {
        onEachFeature: pop_Alle,
        pointToLayer: doPointToLayerAlle
        });

    var searchName = new L.Control.Search({layer: L.featureGroup([json_AlleJSON]) ,
textPlaceholder: 'Bitte einen Suchbegriff eingeben',
propertyName: 'NAME',
circleLocation: false,
markerLocation:false,
zoom: null});

var marker;
searchName.on('search_locationfound', function(e) {

e.layer.setStyle({radius: 5.0, fillColor: 'red', color: 'red', fillOpacity: 1});

   marker = new L.Marker(e.latlng);
   marker.on('click', markerOnClick);
   map.addLayer(marker);
});

map.addControl(searchName);

    function markerOnClick(e) {
    e.layer.openPopup().openOn(map);
    };

How can I open the popup of the feature by clicking its marker?

I thought it might work with e.layer.openPopup().openOn(map); but unfortunately it doesn't

Edit:
Tried both and its not working:

marker = new L.Marker(e.latlng).bindPopup(popupContent);

marker = new L.Marker(e.latlng).setContent(popupContent).open(onMap);


texts: Performance of google map tile creation processes

texts:
I know that question is rather vague but please bear w/me.  I'm trying to get an idea of what sort of product performance - specifically timing - people have seen for various methodologies they have used to create google/bing map tiles.  There are a slew of methods for how to do this (e.g. gdal2tiles, FME, maptiler, etc).  An initial attempt at simply taking a large PNG and creating tiles using imagemagick, on a pretty decent linux server, yielded some pretty long processing times and so I wanted to see what other people are using in production.  New tiles would need to be generated at least daily and so turnaround time on this is pretty critical.

The only real requirement is that it can run on a linux server.  Obviously, free is better but I don't want to restrict myself to that.  The input can be raw gridded/raster data or a large image.  The output needs to be image tiles capable of being used as-is in google or bing maps.

Just for the sake of comparison, I'll say that the timings should be for google map's zoom level 7.

I appreciate everyone's help and again I want to apologize for how vague this question probably seems.

UPDATE: As far as the inputs, I currently have multiple (raw) data sources in various formats: netCDF, GRIB, GRIB2.  In addition to the raw data itself, I also have the ability to generate really large images of that data which could then be sliced/tiled.

Ideally, I would just be chopping the image up but I am willing to try whatever will get me the fastest results.


texts:
I want to to show an image when the mouse hovers on a marker or polygon, and on click show a popup with property attributes information from GeoJSON.

I can manage to show attributes on click but I am not able to show image on hover.
However, the image can be viewed in console.log but it's not appearing on the screen.

Below is what I want it to look like.


texts: Exporting a polyline as KML from AutoCAD and already renamed to its corresponding attribute name

texts: I have multiple closed polylines in AutoCAD with corresponding lot names. How would I go about exporting these polylines and convert them to KML with their names edited to their correct lot names? For example, I have a closed rectangular polygon in AutoCAD that is labeled as Lot 1 (what I meant with labelled is that I wrote a text as "Lot 1" inside that polygon). I like that when I open it as KML in Google Earth, it is not named as "Untitled Polygon" but as "Lot 1". So how do I do this in autocad and in QGIS?


texts: How to apply the orient function on the geometry of a GeoPandas dataframe

texts:
Given a GeoDataFrame where the geometry column contains both Polygons and MultiPolygons.

How can I force the orientation (left-hand-rule or right-hand-rule) on these geometries?


texts:
I was wondering if there is a plugin, or a field calculation in QGIS to calculate a points rotation angle.  ETGeoWizards has a tool that allows you to calculate a point layer rotation angle from a source line layer.  Does anyone know of a similar process within QGIS?

Thanks...


texts: Determining area size of grid (reference grid) in ArcMap

texts:
How do I determine the area size (square kilometres) of a customised reference grid in ArcMap 10.7?
I have tried Data Driven Pages and it is not working for me. The idea is to know the area size of each grid  I made.


texts: Creating raster where each cell records distance to sea?

texts:
I want to create a raster with a 25 metre × 25 metre resolution, where each cell contains the distance to the nearest coastline, as calculated from the center of the cell. To do this, all I have is a shapefile of the coastlines of New Zealand.

I have tried following Dominic Roye's tutorial for doing it in R which works... kind of. It is fine down to about a 1 km × 1 km resolution but if I try to go any higher the RAM it requires well exceeds that available on my PC (~ 70 gb of RAM required) or any other that I have access too. In saying that, I think this is a limitation of R and I suspect that QGIS might have a more computationally efficient way of creating this raster, but I am new to it and I can't quite figure out how to do it.

I have tried following Creating raster with distance to feature using QGIS? to create it in QGIS but it returns this error:

  _core.QgsProcessingException: Could not load source layer for INPUT: C:/..../Coastline/nz-coastlines-and-islands-polygons-topo-150k.shp not
  found

and I am not sure why.

Does anyone have any suggestions of what might be going wrong or an alternative way of doing this?

Edit:

The raster I am hoping to produce would have about 59684 rows and 40827 columns so that it overlaps with the annual water deficit raster from LINZ. If the raster that is produced is larger than the annual water deficit raster I can snip it in R though...

One thing that I think might be a potential issue is that the shapefile of NZ's coastline has a large amount of sea between the islands, and I am not interested in calculating the distance to coast for these cells. I really only want to calculate the values for cells that are include some slice of land. I am not sure how to do this, or if it is an actually problem though.


texts:
After reading Choropleth Mapping Using Quantile Method in CartoDB I decided that it was time for me to understand CartoDB's quantification methods. As I understand it, quantiles should contain an equal number of rows, equal intervals should place the cutoffs at ... equal intervals, and jenks should do some math and come up with groups that have the smallest standard deviation, more or less. I'm skipping heads/tails for now.

I created a slightly goofball map with 18 points, each with an age value from this list:

[11, 11, 12, 13, 14, 16, 19, 20, 20, 20, 23, 25, 26, 26, 31, 37, 38, 38]

I mapped them and then used the choropleth wizard to look at the implication of the various quantification methods. I expect the Quantile method to break the list into groups of six, which it does:

[11, 11, 12, 13, 14, 16]
[19, 20, 20, 20, 23, 25]
[26, 26, 31, 37, 38, 38]

I do not expect to see the logic in the Jenks method at a glance, but indeed, the Jenks groups each have a smaller standard deviation than the quantile groupings.

[11, 11, 12, 13, 14, 16, 19] SD = 2.928
[20, 20, 20, 23, 25, 26, 26] SD = 2.854
[31, 37, 38, 38] SD = 3.367

The SD for the quantile groups is 1.9, 2.3, 5.8 (At least per Calc's =STDEV function.)

So far so good. When I select Equal Intervals I was expecting cutoffs 8 or 9 years apart (so 19 and 28 or 20 and 30 but the CSS places the cutoffs exactly where the quantile cutoffs were, giving groups that are 5, 8, and 12 years long -- not very equal:

#interval_new [ age <= 38] {
   marker-fill: #F03B20;
}
#interval_new [ age <= 25] {
   marker-fill: #FEB24C;
}
#interval_new [ age <= 16] {
   marker-fill: #FFEDA0;
}

Am I misunderstanding how equal intervals work? Or is this a bug?


texts: Writing conditional python script to delete records?

texts:
I am new to python programming. I am wondering if someone could tell me how to write a python script that would find all records in an attribute table that have identical values in one field, then from this subset find the records that are identical based on a second field, then from this subset of records delete record(s) based on conditions set for a third field. So for example, if five records have the same values in field "C", and two of these records have the same values in field "A", but have different values in field "B", I want to delete just one of these records, specifically the one that has a "null" value in field "B".

The attribute table has ~5000 records, I need to delete about ~400 of these records b/c they are identical to other records in this table.

In the example shown below, I would want to use Python to find the records with identical values in "FEAT_SEQ*" field, then from this subset find records with identical values in the "SHAPE_len" field, then from this subset delete the record with a "Null" value in the "edited" field. In this example, I would keep three records, and delete the record that is second from the top ("IN_FID*" = 23)

I am working in ArcGIS 10.1


texts: How to determine if a LYR file has more than one Layer using ArcObjects?

texts: In using arcpy, I have found that there are just a few lines of code needed to get a list of layers in a LYR file if more than one layer exist.  How can I accomplish this same feat using ArcObjects in C#?  I know that ILayerFile is the way to start this off and I can open a specific layer file using ILayerFile.Open() but how can I get back a list of layers in this layer file if more than one exist?


texts: Labeling only maximum value when zoomed out or print in smaller scale

texts:
I would like to prioritise the maximum values to be shown on the map. In other words, when the label collides it should drop out the minimum values.

In the shown examples, I would like to have priority for -5.86 over -5.91 when I print/zoom in a smaller scale.


texts:
I've the floor plan for a building in PDF. I'm trying to upload this into Mapbox as a tileset but I'm having a lot of confusion regarding the conversion.

How should I approach this?

P.S. since I'm new to GIS, please provide detailed explanations


texts: How to create dissolved buffers using QGIS virtual layers?

texts:
To create dissolved buffers (single geometries) in PostGIS I use ST_Dump and ST_Union:

SELECT
    row_number() over() AS gid,
    sbqry.*
FROM (
    SELECT
        (ST_Dump(ST_Union(ST_Buffer(input_layer.geom, 100)))).geom::geometry(Polygon, 31256) AS geom
    FROM input_layer
) AS sbqry;

As far as I know QGIS virtual layers are based on SpatiaLite. So at least ST_Dump isn't supported. Is there a SpatiaLite equivalent to the query posted above?


texts: Why does QGIS add decimals when a csv is saved as a shapefile?


texts: Is there anyone who can help me with SAR data. I am a beginner with this. Any tutorial/links will help.


texts: I want to to find the social welfare and the population density for Europe and Africa with night-lights. Is there any site where I can find images with good resolution (I only think VIIRS)??? I have never worked with night-lights so I don' t know much.


texts: JTS GeometryClipper returns wrong results

texts:
I have envelopes and a polygon and I want to split that polygon using the envelope geometries. My code is like this:

val arr = new ArrayBuffer[Polygon]
for (h in polygons){
    val ex =  h.envelope
    val cellEnvp = new Envelope ( ex.Xmin, ex.Ymin, ex.Xmax, ex.Ymax )
    val clipper = new GeometryClipper ( cellEnvp )
    val clip = clipper.clipSafe( h.jtsGeom, false, 0.001)
    if (clip != null && !clip.isEmpty) {
      arr += clip.asInstanceOf[com.vividsolutions.jts.geom.Polygon]
    }
  }

The thing is that the result is OK, but it adds in each polygon a little part from the envelope. Why...?

UPADATE: if I use true, I see the above results
 1. Now, it returns polygons and also Multipolygons
 2. I loose some clips

Results: The red polygons are the one that missing


texts: Converting Lat/Lon to Pixel-XY on cropped map?

texts:
I have a GPX-Track containing LatLon coordinates based on WGS84. The track should be drawn on a cropped map, and for this, I need the pixel coordinates. The map tiles are from OSM

I'm not quite sure how to do this, but I figured out a way (without implementing it yet, since I don't have the formulars right now).

Convert LatLon to Spherical Mercator which is used by OSM
Linear transformation spherical mercator -> pixel coordinates.
For step 2., 2 identical points are needed. I guess I can use the map's corner points.

Is there a better way to do this?


texts:
I'm working on a map showing on which days garbage is collected in different areas. That means it's pretty tempting to colour code each zone with a different colour representing a day of the week.

Are there any conventions to follow? Is 7 colours too many?

Here's what it looks like at the moment:

My colours so far follow the rainbow from Monday to Friday (red, yellow, green, blue, purple) then light grey and dark gray for Saturday and Sunday.

Here are the numerical values in the image, for interest. (CartoCSS, Tilemill.)

@s: 60%;
@l: 70%;
@monday:    hsl(  0, @s, @l);
@tuesday:   hsl( 60, @s, @l);
@wednesday: hsl(110, @s, @l);
@thursday:  hsl(200, @s, @l);
@friday:    hsl(250, @s, 80%);
@saturday:  hsl(  0, 0,  70%);
@sunday:    hsl(  0, 0,  70%);

So the weekdays are all the same luminosity and saturation except Friday, which I found seemed too dark, for some reason. I originally spaced the hues equally around the spectrum, but found that there were two blues too similar, and several of the colours just weren't attractive.


texts: Get a point value with psycopg2 and PPyGIS

texts:
I'm trying to setup a basic working postgis setup with python ppygis package.

>>> import psycopg2
>>> import ppygis
>>> connection = psycopg2.connect(database='spre', user='postgres')
>>> cursor = connection.cursor()
>>> cursor.execute('CREATE TABLE test (geometry GEOMETRY)')
>>> cursor.execute('INSERT INTO test VALUES(%s)', (ppygis.Point(1.0, 2.0),))
>>> cursor.execute('SELECT * from test')
>>> point = cursor.fetchone()[0]
>>> print point
0101000000000000000000F03F0000000000000040
>>>

I should have got a python object with separate X and Y coordinate. Something like

>>> Point(X: 1.0, Y: 2.0)

What am I doing wrong? I'm following http://www.fabianowski.eu/projects/ppygis/usage.html#installation

For posterity. I couldn't make it work, so to get point data this was what I did,

>>> cursor.execute('SELECT ST_AsGeoJSON(geometry) FROM test')
>>> point = cursor.fetchone()[0]
>>> print json.loads(point)['coordinates']
[1,2]


texts: How do I delete lines from the QGIS canvas which no longer are attached to a shapefile?

texts: Using QGIS 2.18.1, I started to create a shapefile layer with a line feature, but then decided not to proceed.  I first tried to select and delete the lines but couldn't seem to select them with the select tool.  So I then deleted the shapefiles using File Explorer which seemed to go as expected.  However, the lines are still on the canvas, and appear to be impossible to delete, and of course are no longer associated with a layer .  How can I clear the canvas of these lines?


texts: display issue in qgis 2.14 : features appear and disappear

texts:
I have several virtual layers in my project. One of them contains 6 lines. Depending of the scale, at least one line disappears from the canvas for an unknown reason.
I heard that sometimes, creating spatial index solve the problem.
I tried to create an index for this virtual layer according to spatialite sql (select CreateSpatialIndex('table', 'geometry')) but it fails (returns 0), i don't know the reason why.
Question : why does the index fail and/but is the display problem really linked to the spatial index ?
Here are two pictures with two differents scales.

edit : added shared data
https://drive.google.com/file/d/0B1yuAbW2T79iU0hDUFpFanZWLTg/view?usp=sharing

SQL :

select t1.* from "Tours_cupidon_unobs_ordre_rendu"
as t0,"Tours_cupidon_unobs_ordre_rendu" as t1
where st_distance(st_centroid(t0.geometry), st_centroid(t1.geometry))
<4 and t0.num = 314 AND t1.typFait = 'sep'

select t1.num, make_line(centroid( t0.geometry), centroid(t1.geometry))
as geometry from "Tours_cupidon_unobs_ordre_rendu" as t0
, "Tours_cupidon_unobs_ordre_rendu" as t1
where t0.num = 314 and t1.num in (324,325,334,347,348,349,3142,3153,3162)AND t1.typFait = 'sep' and t1.typunobs = 'fait'


texts: Leaflet - Check if polyline intersects polygon

texts:
What is the best way to determine if a polyline (not just the vertices) intersect a polygon? Or checking if the polygon contains any of the polyline?

I am checking within a function where my polyline is an array of points like below.

var myPoints=[[42.414649077488, -71.13123893737793], [42.413286665076754, -71.13115310668945], [42.41220938779894, -71.13089561462402]]

If necessary I'd be happy to convert to a polyline feature like below before checking

var myPolyline = new L.Polyline(myPoints,{
        color: '#38c3de',
        opacity: 0.9,
        weight: 10,
      })


texts: Calculating Euclidean distance in QGIS 2.18.27

texts: I want to calculate Euclidean distance of animals trails from roads. Is this possible as both layers or vectors? Or is there another way to calculate distance my animal trail network is from a road?


texts: Calculating Datum Corrections to Overlay pre-NAD27 USGS Maps?

texts:
I am working with some older (ca. 1900) USGS Topographic Maps. These Raster maps are georeferenced  with a CRS of NAD27 and poly projection (polyconic). However, they are offset and require easting/northing corrections in the georeferenced source data to overlay on modern maps.

On some maps, in the 1920's, Someone had calculated these offsets for the original editions, perhaps to celebrate the creation of NAD27, who knows.

I am looking for some pointers on how that was done.

I want to know if there is a way to calculate this correction generally/automatically. , or alternatively if there is a standard set of datums to use to account for the easting/northing correction in these old maps.

Example

San Francisco 1916 15' x 15', Warning >10MB

This source data is already referenced to NAD27, the CRS as encoded in the file is usually interpreted by QGIS as something like

+proj=poly +lat_0=37.75 +lon_0=-122.375 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs

While the projection is good, it is not perfectly correct, when used as-is the data is quite offset

Luckily, on some editions the following can be found

Which I Interpret as Easting and Northing corrections to the CRS definition.

+proj=poly +lat_0=37.75 +lon_0=-122.375 +x_0=-97.536 +y_0=-201.168 +datum=NAD27 +units=m +no_defs

The result is quite good.

Additionally, I have also encountered the following note

By extracting the measured projected distance (in QGIS) between the solid frame and the dashed frame

I can extract x_0 and y_0 and use it like above

Sometimes, however nothing is provided and I calculate this offset by hand by doing some projected distance measurements of good reference points.

In all cases its a labor intensive process.

Last Word

Someone has calculated these offsets before, so I am looking for some pointers on how that was done.


texts: Is there a standardized or "most used" dummy Z value?

texts:
Creating and importing both 2D and 3D data, I have many times encountered the situation where I have no Z value for a coordinate set, that the value of a Z coordinate seems out of range (like -99, -9999, -inf or similar) or that I need to create a dummy Z coordinate.

I know that the answer to my question is:

  "Just use a value which is definitely
  out of range in your case."

But that answer set aside I wonder if the GIS community has a standardized or most frequently used value for a dummy Z coordinate?


texts: Pansharpening problem in ArcGIS: failure when using TOA calibrated Landsat

texts:
I'm hoping to use the "Create pansharpen raster dataset" tool in ArcGIS to pansharpen Landsat 7/8 images. This tool works nicely and produces the expected pansharpened image when the input is in DN form. However, when I input TOA reflectance values (both multispectral and pan) instead of DN's, the function fails, and produces an all 0 tiff.

Does anyone know why that would happen?


texts: I want to perform JOIN TOOL for all DBFs generated from ZONAL STATISTICS AS TABLE into one table. However, it seems that DBF has a maximum of 255 fields while my dataset has more than 500 fields/columns. So, my question is how can I divide/separate the INPUT FOLDER such that ZONAL STATISTICS will perform 250 tables at a time? All else is working other than this issue. I am using ArcGIS 10. Thanks.


texts: Overlay image in Mapbox gl doesn't work in Chrome

texts:
I'm trying to prepare image overlay layer in mapbox gl with city plan as background. It works perfect in Firefox but in Chrome I see only dark grey polygon. I use 4Mb .png file.

Firefox Desktop - ok

Firefox Mobile- not ok

Chrome Deskotop -not ok

Firefox Mobile -not ok

       map.addLayer({
                "id": "overlay",
                "source": {
                    "type": "image",
                    "url": "7.png",
                    "coordinates": [
                        [30.0112567, 51.4361221],
                        [30.187844, 51.4361221],
                        [30.187844, 51.3696605],
                        [30.0112567, 51.3696605]
                    ]
                },
                "type": "raster",
                "paint": {
                    "raster-opacity": 0.85
                }
            });


texts: Getting the reflectance images in red band and near-infrared band respectively corresponding with the max NDVI

texts:
How to get the reflectance images in red band and near-infrared band respectively corresponding with the max NDVI that I calculate from the reflectance image collection in their red bands and near-infrared bands? The following code shows how I get the image of max ndvi.

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T2_SR')
    .filterDate('2016-05-01', '2016-9-30')
    .filterBounds(fc);   // fc is the polygon of study area.
// A function to compute NDVI.
var NDVI = function(image) {
  return image.expression('float(b("B4") - b("B3")) / (b("B4") + b("B3"))');
};   //B4 is the reflectance in the near-infared band, and B3 is the reflectance in the red band.

// Reduce the collection by taking the max.
var max_ndvi = collection.map(NDVI).max();


texts: Appending polygons to existing feature class in ModelBuilder?

texts:
I have a shp file with lots of polygons. I would like to create a buffer zone around the first polygon and intersect it (the buffer) with an other feature layer, than store the result in an existing feature class and iterate this process with Iterate Feature Selection (by ID) till the last polygon. The model runs without any error, I see that it creates every buffer, but the final layer contains only the last one.

I would like to use the layer of buffer zones just as a temporary layer and overwrite it in every iteration and just store one feature layer with the intersections of buffers. How can I append the intersections of buffer zones to the final layer without overwrite the previous result?


texts: How to pass expression to where clause in ArcPy?

texts:
I want to pass a user defined variable into a where clause

Syntax
SelectLayerByAttribute_management (in_layer_or_view, {selection_type}, {where_clause})

#user input
MTRsel = arcpy.GetParameterAsText(0)

#put user input into where clause
arcpy.SelectLayerByAttribute_management("Townships","NEW_SELECTION", "MTR" = + MTRsel)

I have tried many variations of this but it will not accept it, what am I missing


texts: Extract values from Grib-files at several specific locations

texts:
I have got a GRIB file from ECMWF for a small area in the Antarctic. It contains hourly data for one month of multiple variables (eg. wind, seaicecover, sst etc.). I want to extract this data to .txt files (every file for a different variable) at several specific locations...meaning that i have around 1000 sampling locations (stored in a text file with rows for Longitude, Latitude and ID) , for which i want to get the environmental data of the GRIB. Is there a way to get this done without doing manually? In the end a file (eg. windspeed) should have all the windspeed-values for given locations at the different times. I have already tried to use cdo and wgrib2 for example but neither of them provided me with a fitting solution or i didn't found the right tools. I am using Windows 10 (64 bit) and i've got Python 3.8 and Ubuntu 20.04 LTS installed. If you would need more specific information please let me know. I cannot include any code here because i have absolutely no idea how to approach this problem.
But nevertheless here is an example of the Format of the sampling locations:
Longitude  Latitude fid
1 -50.87164 -61.57444  R1
2 -60.65056 -61.30030  R2
3 -50.46169 -61.01206  R3

Do you have any idea which tool i could use for this?


texts: EarthEngine - Layer error: Image.parseExpression: Expression parse error

texts:
I am trying to create a layer using conditional statement in Earth Engine. I have created a stack of two layers as follows:

var NDVI_2000_2001_2002 = NDVI_2001_2000.addBands(NDVI_2002_2000);

Now, I want to create a layer with three classes, where common pixels having values greater than 0.2 in both the layers (NDVI_2001_2000 and NDVI_2002_2000) will be assigned a new value 3; values less than -0.2  a new value of 1 and everything else will be assigned a value of 2. I am using the conditional statement as follows:

var Change_2000_2001 = NDVI_2000_2001_2002.expression(
    "(b('NDVI_2001_2000', 'NDVI_2002_2000')) > 0.2 ? 3" +
      ": (b('NDVI_2001_2000', 'NDVI_2002_2000')) > -0.2) ? 2" +
        ": 1"
    );

I am getting the following error:

  Layer 2: Layer error: Image.parseExpression: Expression parse error at
  character 97: '(b('NDVI_2001_2000', 'NDVI_2002_2000')) > 0.2 ? 3:
  (b('NDVI_2001_2000', 'NDVI_2002_2000')) > -0.2) ? 2: 1'
                                                                                ^.

I have not been able to identify where exactly the correction needs to be done.  Kindly advice. Also, any alternate ideas to create the layer is also welcome.


texts: Should I use Basemaps or Layers to overlay two maps?

texts:
I'm trying to create a map that essentially has two maps on top of each other (chosen from many map choices) at the same time. One of maps will have an opacity slider. Should I be trying to use two basemaps? or should I be adding them as layers? What is the best way to go about doing this?

It should let there be an effect like this with the first map showing through the second.


texts: Table Join on distance measures returns NULL

texts:
I followed every step of this tutorial but when I join the distToRoads and distance in m I get a table distToRoads which contains just the null vaues, why is that?

Tutorial:  https://www.qgistutorials.com/en/docs/performing_table_joins.html


texts: Remove .sr lock on shapefile after removing layer from map using ArcObjects and C#

texts:
I am having difficulty figuring out how to remove .sr lock files from shapefiles using ArcObjects and C# after removing the associated layer from the map.  I notice that using the straight UI in ArcMap that removing a layer will remove the lock file, but I cannot get this to happen programmatically.  I have tried everything that I can think of, including calling Marshal.FinalReleaseComObject on all variables.  The following is some paired down test code that simply loads a layer into the map and then immediately removes it.  This should remove the lock files, but it does not.  Am I missing something?
// Loop through the selected shapefiles and load into the map
pGxObjects.Reset();

while ((pGxDataset = (IGxDataset)pGxObjects.Next()) != null)
{
    IFeatureLayer pFeatLayer = new FeatureLayer();
    pFeatLayer.FeatureClass = (IFeatureClass)pGxDataset.Dataset;
    pFeatLayer.Name = pFeatLayer.FeatureClass.AliasName;
    pMxDoc.FocusMap.AddLayer(pFeatLayer);

    // Do something with the layer

    pMxDoc.FocusMap.DeleteLayer(pFeatLayer);
    Marshal.FinalReleaseComObject(pFeatLayer);
    Marshal.FinalReleaseComObject(pGxDataset);
}

Marshal.FinalReleaseComObject(pGxObjects);
GC.Collect();
GC.WaitForPendingFinalizers();


texts: Where to start with great circle maps in CartoDB?

texts:
I've pulled together a series of maps in R based on Nathan Yau's excellent Great Circles tutorial and a nice zipcode package The maps show how many people from many zip codes descended on a single zip code. Here's one example:

I'd like to reproduce the maps in CartoDB (which uses PostGIS, I believe), but I'm not sure where to start.

My read of what I was doing in R was generating a matrix called "inter" that was full of lat/long points, drawing the line connecting those points and then starting over. So one solution that occurs to me is to have R store each line as a map line and take that series of lines to CartoDB. But I'm wondering if there's a reasonable way to just generate the maps in CartoDB directly, given an end point and a series of starting points.


texts: How to do a search in a line to a specific point

texts: I'm using QGIS 2.0.1 Dufour. I'm working with a shape file made by a single line that represents a road with 30km. In that line I want to attach pdf files containing information about some points in that road. I would like to be able to do a search like: 'km = 2,300' and find the files with the information about that point. Is that possible?


texts: Join on range query using Esri's GIS Tools for Hadoop?

texts:
I am getting myself acquainted with Esri's GIS Tools for Hadoop. Working with Hive.

One of my tables, A, has user lat-lons. Another, B, has landmarks and their lat-lons. How do I join two tables so that I get all landmarks names and their locations from B within a given distance (radius) from each location in A? How to replicate a k-NN query, if possible?

Can I do it with Esri's GIS Tools for Hadoop?

I have heard of Hadoop GIS and SpatialHadoop, but there are limitations - cannot get them working.


texts:
I am trying to use the Raster Calculator to multiply 2 rasters together.  They were created by using the raster calculator to keep values exceeding a threshold in two other rasters.  I have tried numerous combinations but always end up with zeros everywhere.

I have 2 rasters A1 and B1 of the same area.  I want to locate areas where A1 > 5 and 30 <= B1 <= 45.

I could not figure out how to use AND or OR commands in the Raster Calculator (the RasterCalc plugin suggested elsewhere doesn't appear to still be available).

I created two new rasters A2 and B2, using the commands

("A1" > 5)    {saved as A2}

("B1" <=45 )* "B1"   {saved as B2a}
(B2a >= 35)          {saved as B2}

I now have 2 rasters that are 0,1 values (or rather 0.999000).  I want to multiply these together.  I have tried wrapping them in float(), or using float with the above commands to see if that made a difference (it didn't).

If I look at the rasters/use the inspection tool  I can see that there are regions where both A2 and B2 are 1 - yet any multiplication gives only 0 as the output, any addition gives only 0 or 1.

I have not been able to find any reason as to why this might be happening.

If I simply multiply the two original rasters A1 and B1 this works as expected, so it does not appear to be multiplication that is the problem.


texts: Dissolving in ArcMap gives "Invalid Topology [Duplicate Segment.]"

texts:
I am attempting to use the dissolve tool on a shapefile, but it fails with the error:

ERROR 999999: Error executing function.
Invalid Topology [Duplicate segment.]

I tried checking and repairing geometry until no more geometry errors exist. I also tried the fix topology error tool, but the topology edit tool says that no topology exists. Any ideas why this might be happening?


texts:
I have made the UI of a qgis plugin using qt creator. What I want is that the plugin window should stay on top of the qgis main window only but not other other applications. I have got this code:

    QtGui.QDialog.__init__(self, None, Qt.WindowStaysOnTopHint)

But the problem is that after implementing this code the plugin window stays on top of every other application as well.


texts: Creating table in ArcToolbox Script with ArcPy

texts:
I want to create a Python Script for ArcGIS, using arcpy.
The script would have several user input fields. I would like to create a table based on one or two of the user input (the path and name).
I am trying :

import arcpy
import os
from arcpy import env
from arcpy.management import CreateFeatureclass
workspace = arcpy.GetParameterAsText(0)
poly = arcpy.GetParameterAsText(1)
poly_nou = arcpy.GetParameterAsText(2)
input_areas = arcpy.GetParameterAsText(3)
direction = arcpy.GetParameterAsText(4)

env.workspace = workspace
out_path = os.path.split(poly_nou)[0]
name =  os.path.split(poly_nou)[1]

arcpy.CreateTable(out_path, name)

I have set up the Script in ArcMap to look like this:
1. Workspace type
2. Feature class
3. Table
4 and 5 Strings

I want the script to read 3 (possibly also 1) and create a table.

I get the error:

  Traceback (most recent call last):
    File "C:\Users\user\Documents\Proiecte\Python Toolbox\Divide_percentage.py", line 16, in
      arcpy.CreateTable(out_path, name)
  AttributeError: 'module' object has no attribute 'CreateTable' Failed to execute (Script1).


texts: Which web application framework is easiest for building website displaying lot of geospatial data?

texts:
I am building a simple website that displays a lot of geospatial data (which is routed through Geoserver). I was trying to use Drupal-6 to build the website component (non geospatial descriptions) but had plenty of issues dealing with openlayers and other modules as well as deploying it properly.

My goals are simple - to display geospatial data (with Google maps background) and a normal website with a big menu and different explanatory content.

Which web application framework is easiest (for a non web programmer) to achieve these two goals based on your experience


texts: Legend Key for discrete values in spplot

texts:
I have a shapefile with various attributes. I want to plot the shapefile using `spplot'.

I first read the shapefile

A <- readOGR('myshapefile.shp')

I am using the following command:

p <- spplot(A, c('lulc', 'lulc_2', 'lulc_3',
             'lulc_4', 'lulc_5', 'lulc_6', 'lulc_7'),
        col.regions=myColors,
        col='transparent',
        names.attr = c(1998, 1999, 2002, 2003, 2009, 2014, 2017))

Each of the fields / attributes contains 4 values, say 1, 2, 3, 4 and correspondingly need to be labelled A, B, C, D in the legend, along with the color key. How to achieve this? For example, I have put the image of the legend only.

Small update:

I tried to modify the spplot command to:

p <- spplot(A, c('lulc', 'lulc_2', 'lulc_3',
             'lulc_4', 'lulc_5', 'lulc_6', 'lulc_7'),
        col.regions=myColors,
        col='transparent',
        names.attr = c(1998, 1999, 2002, 2003, 2009, 2014, 2017),
        as.table = T,
        colorkey = list(
          right = list( # see ?levelplot in package trellis, argument colorkey:
            fun = draw.colorkey,
            args = list(
              key = list(
                at = c(1,2,3,4), # colour breaks
                col = myColors, # colours
                labels = list(
                  at = c(1,2,3,4),
                  labels = c('A', 'B', 'C', 'D')
                )
              )
            )
          )
        ))

Even then this is not working.

(see https://www.rdocumentation.org/packages/sp/versions/1.3-1/topics/spplot)


texts:
I have a shapefile of a city. My problem is that in the border, cells of the grids are not a complete square (500*500 m), while I need them to be complete as I need their centroids too. As it is shown for example in the picture, all cell in the border should be complete( not limited to the city border).
The steps I followed is that through fishnet I created grids and then clipped it to the city border.
How I can fix this issue?


texts:
How would one to check to see what was dissolved using the Dissolve tool in Arcmap 10.2?

Other than manually checking to see what polygons were dissolved by toggling back and forth between the dissolve input and output layers is there a method, tool, or way to visually verify what was dissolved?


texts: How does one install QGIS from previously downloaded OSGeo packages

texts:
I have downloaded the OSGeo Installer (on Windows 7 x64) and run it, selecting the Desktop Express Install option.  All of the files were downloaded and the installer reported a successful install, however whilst I could find the downloaded files, I could not find the installed instance.  An OSGeo4W directory was created but contained only a few MB of files.  I have tried to rerun the installer and it only confirmed that the previous installation was successful and that no files needed downloading.  All the downloaded files were still in tarballs so it appears that the installation process after downloading had not taken place.  QGIS has previously been on my system but all previous and development versions were uninstalled using the provided uninstallers.

How to I bring about an installation process that installs desktop QGIS using the files downloaded by the OSGeo Installer?


texts: Generate a list of intersected layers QGIS

texts:
Apologies, im very new to QGIS and there's a chance what I'm asking for isn't actually what I require. Let me explain.
I'm working with a project with multiple layers. I have a layer (let's call it layer 1) with multiple polygons outlining site boundaries. What I would like to know is if it is possible for QGIS to generate a list of all the layers that a specific polygon in layer 1 intersects with?
I know how to use a spatial query but all this does is highlight a load of polygons that intersect with a particular layer.


texts: Is there a tool that allows the same functionality of extend/trim feature that was available in 9.3?

texts:
In ArcGIS 9.3 there was a tool called Extend/Trim Features that allowed you to select multiple segments and extend or trim them to a feature.

In 10.1 that seems to have been moved to the Production Mapping Extension at a hefty price. It was replaced by two tools that don't seem to have the same functionality as the previous versions tool.

My question: Is anybody aware of an add in or way that I could create this functionality with python or a custom tool?

EDIT

The old tool can be found referenced HERE, there is also a screenshot of the tool location from 9.3:

Essentially, it did both of the functions that are shown in @Jakub's answer.

It looks like the same tool from 9.3 was moved to the Production Mapping toolbar in 10.1.

There seems to be an ESRI Idea out there suggesting they bring it back, and that I'm not, in fact, imagining that it used to be this way.


texts: How to get the job status of gis.item.export or make it execute synchronously?

texts:
I'm using the ArcGIS API for Python to export a feature service to a FGDB and then download it.
My problem is that the export function appears to be asynchronous even though I set wait=True.
The returned item from gis.item.export() has a jobid so this shouldn't be a problem if I could poll the service for job status but I don't see a URL to do this.
Without a synchronous export or a job status I am user time.sleep(60) to slow things down but that not exactly ideal
So how do I get my export to execute synchronously, or how do I query the job status?
def Get_AGOL_Data_All2(itemID):
    """
    PARAMETERS:
      itemID (str) = Ther Portal itemID of the feature service to download.

    RETURNS:
      None

    FUNCTION:
      Uses the arcgis api for python library to connect to Portal and copy a feature service to a
      FGDB. Portal returns a zipped file which then has to be unzipped and the contents moved to the staging location

    NOTE:

    """

    print ('Starting Get_AGOL_Data_All()')
    # #anon_gis = GIS()
    mygis = GIS(username="userid",password="****")

    service = mygis.content.get(itemID)
    itemDesc = service.export(title="DamageReports", export_format= 'File Geodatabase', parameters=None, wait='True')
    time.sleep(60)
    fgdb = mygis.content.get(itemDesc['exportItemId'])

    fgdb.download(r"D:\temp")

    with zipfile.ZipFile(os.path.join( r"D:\temp" , fgdb.name),'r') as zip_tar:
        zip_tar.extractall(r"D:\temp")
    print ('Finished Get_AGOL_Data_All()')
    arcpy.Rename_management(os.path.join(r"D:\temp", zip_tar.filelist[0].filename.split("/")[0]),"DamageReports.gdb")

    return


texts: Legal Requirements for using Drones/UAS/UAV for GIS

texts:
Can anyone give me a better idea of the specific requirements for commercial use of drones for GIS work? What would be required of the pilot?

Edit: I'm also looking for more insight into when the Small UAS Notice of Proposed Rulemaking goes into effect. I'd interested in operating smaller (<55 lbs) drones which would hopefully not require an FAA license.


texts: Cannot drag polygon features in ArcMap 10.1

texts: In ArcMap 10.1, I can 'Start Editing' a polygon layer and select a feature with the 'Edit tool'. If I hold down the left mouse button, I get the four way arrow cursor but then I am unable to drag the polygon anywhere. I can however edit vertices and change the relative position using the 'Move" command. Same issue with all shapefiles and projects.  This functionality works fine on the different machine with ArcMap 10.1. I assume that something in the settings is preventing this type of editing.


texts: Understanding the polygon handler function in Openlayers

texts:
Using openlayers, I want to be able to draw one polygon on a map. After one polygon is drawn, I want an event to fire that stops any more from being drawn until a button is pressed. My current code allows the user to draw more than one polygon and I'm trying to figure out how to limit it. This is my current code:

// Add drawing layer
analysisLayer = new OpenLayers.Layer.Vector("analysisLayer");
map.addLayer(analysisLayer);

// add polygon drawing control
polygon = new OpenLayers.Control.DrawFeature(analysisLayer, OpenLayers.Handler.Polygon);
map.addControl(polygon);

What I tried doing was adding a callback to the Control.DrawFeature function as follows:

// add polygon drawing control
polygon = new OpenLayers.Control.DrawFeature(
    analysisLayer,
    OpenLayers.Handler.Polygon, {callbacks: { "finalize" : function(){alert('poly made')} } } );
    map.addControl(polygon);

This does not work, but it does not give me an error. I'm not sure how exactly to structure this function. I did find the function finalize in the documentation which is what I thought I needed, but I'm also not sure how to structure the callback.

EDIT:
I'm starting to understand how these callbacks work. For example this detects that a feature has been added (using @Vadim's advice):

// detect polygon events
function newPolygonAdded (evt) {
        alert('Polygon completed');
        polygon.deactivate(); //stops the drawing
    }

// add polygon drawing control
polygon = new OpenLayers.Control.DrawFeature(analysisLayer, OpenLayers.Handler.Polygon,
                            {eventListeners:{"featureadded": newPolygonAdded}});
map.addControl(polygon);

What I still don't understand is how to limit the drawing layer from allowing more polygons to be drawn.


texts:
I am trying to make a web application that renders 3D buildings (LoD1). I have a PostGIS (2.2.1dev) database which contains the 2D geometry of the buildings and the height. I used st_extrude to make a 3D geometry out of it.

The problem is that I end up with a PolyhedralSurface, while I need a MultiPolygon. I used st_dump to split up my polyhedralSurface into its faces (floor, roof, walls). I end up with a table that has an id-column and a geometry-column.

For each building, I now have multiple rows in the table. The geometry column is of the form:

({1}, 01030000A8...).

The {i} part is the identifier for which face of my building it is, the second part is the actual geometry of this face.

How can I join the multiple rows for one building together and end up with a MultiPolygon?


texts: Create points along a track at a specified straight line distance

texts: I'm using ArcMap 10.5. I have a track that I wish to place 30 points along at 400 meter intervals. The issue I am facing is that these points must be 400m apart measured in a straight line, not the actual distance along the track (see pic). I could do it manually but this would be quite time consuming. Is there a tool or extension that I can use in ArcMap to achieve this?


texts:
I am now want to use subprocess.call with gdalwarp to clip raster from shapefile
(in qgis 1.8.0)

I use raster→projections→wrap(reproject)to generate command_line below

gdalwarp -t_srs "+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0
+ellps=GRS80 +units=m +no_defs" -dstnodata 0 -q -cutline "C:/Program Files
(x86)/Quantum GIS Wroclaw/test/area.shp"
-dstalpha -of GTiff C:/Program Files
(x86)/Quantum GIS Wroclaw/test/wrapM_0.tif "C:/Program Files
(x86)/Quantum GIS Wroclaw/test/wrap"

and following is my code

import gdal
from subprocess import call

crs='"+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0
+ellps=GRS80 +units=m +no_defs"'

mask='"C:/Program Files(x86)/Quantum GIS Wroclaw/test/area.shp"'
input='C:/Program Files(x86)/Quantum GIS Wroclaw/test/wrapM_0.tif'
output='"C:/Program Files(x86)/Quantum GIS Wroclaw/test/wrap"'

call(gdalwarp, -t_srs crs -dstnodata 0 -q -cutline msak -dstalpha -of GTiff
input output)

And Python will give the error below:

  File "<input>", line 1
    call(gdalwarp, -t_srs crs -dstnodata 0 -q -cutline msak -dstalpha -of
GTiff input output)
                            ^
SyntaxError: invalid syntax

How to fix the error here??


texts: ArcGIS Pro: Select by location result varies w.r.t system being used

texts:
We are getting different results if we run the "Select by location" tool in ArcGIS Pro (same version 2.1.1) for same two layers but in different machines. We didn't assign any tolerance level and the units are kept "Unknown". How is it possible ? Does it have to do anything ArcGIS pro settings ?.

Edit: Attached is the screenshot of parameters used. Left the values to default

Does it have to do anything with M tolerance and M Resolution ?


texts: How to show marker on buffer and show the info on click of that marker ArcGIS JavaScript?

texts:
First for creating the buffer the below function is used

function custombuffer(radial,latt,longt,zoomlevel) {

        map.graphics.clear();
        var pt = new Point(Number(longt), Number(latt), new SpatialReference({
                                wkid: 4326
                            }));
                             map.centerAndZoom(pt,zoomlevel);
        var params = new esri.tasks.BufferParameters();
        params.geometries = [pt];
        //buffer in linear units such as meters, km, miles etc.
        params.distances = [radial];

        params.unit = esri.tasks.GeometryService.UNIT_KILOMETER;
        params.outSpatialReference = map.spatialReference;
        gsvc.buffer(params, showBuffer);

    }
    function showBuffer(geometries) {
        var symbol = new esri.symbol.SimpleFillSymbol(
                esri.symbol.SimpleFillSymbol.STYLE_SOLID,
                new esri.symbol.SimpleLineSymbol(
                    esri.symbol.SimpleLineSymbol.STYLE_SOLID,
                    new dojo.Color([0, 255, 0, 0.1]), 2),
                new dojo.Color([0, 0, 255, 0.2]));

        dojo.forEach(geometries, function (geometry) {
            var graphic = new esri.Graphic(geometry, symbol);
            map.graphics.add(graphic);
            //map.centerAndZoom(pt,10);

        });
    }

and then I hit the webservice to get the lattitude and longtitude to
create the marker inside the buffer which is created using above
function.

Here is the code:

  $.ajax({
                type: 'GET',
                url: "app_code/pointdata.php?filtervalue="+filtervalue+"&filtertype="+filtertype+"&latt="+latt+"&longt="+longt+"&radial_dist="+radial_dist,
                crossDomain: true,
                success: function (dat) {
                    var dataPoint = JSON.parse(dat.trim());
                    for (var i = 0; i < dataPoint.length; i++) {
                if (dataPoint[i].longitude && dataPoint[i].latitude) {
                    var pt = new Point(Number(dataPoint[i].longitude), Number(dataPoint[i].latitude), new SpatialReference({ wkid: 4326 }));
                 //   var markersymbol = new PictureMarkerSymbol('images/marker.png', 10, 15);
                   var content = '
';
                    content += "
test	" + ",	test
";

                    var infoTemplate = new InfoTemplate('Point info', content);
                    polygraphlayer.add(new Graphic(pt, markersymbol, attr, infoTemplate));

                 // polygraphlayer.add(new Graphic(pt, markersymbol));
            }
             }

buffer and marker are created successfully but the problem is when I
try to click the marker its simply zoom in the buffer. I want to show
the content info on that marker


texts: Use GetFeatureInfo with OpenLayers and a Bounding Polygon

texts: I want to be able to draw a bounding polygon and have my MapServer WMS layer return every feature within that bounding box. I have it working for a single point, but I would like for the query to return every single feature within a polygon. Is that possible with WMS GetFeatureInfo? I am using OpenLayers in Java with GWT. There is a bridge between Java and the OpenLayers JS.


texts: Removing coincident points with z values greater than minimum

texts:
I am trying to run a spline interpolation on some terrestrial LiDAR data, but it keeps failing (I'm getting error 010235 - returned error from spline process and error 010067 - error executing grid expression).

I think it's due to the fact that there may be coincident points (ie. points with the same xy but different z values).

Can anyone suggest a way to delete the duplicate points with z values greater than the minimum? I can't do it manually because I think there may be a lot of duplicates (I tried the collect events tool, and there were many coincident points).


texts: Cleaning up data in Excel for use in QGIS

texts:
I am trying to clean up a large .txt dataset in Excel.  I think that a macro or another application may be able to clean up the data they way I want but I have limited programming skills.  I need to create fields from different sections in the .txt file so I can then geocode the points.  I am conducting a crime analysis for my local police department and up until now I have been using wildcards and the Find & Replace tool to do this, but I have nearly 2 years worth of police calls.  Fields needed:"Address", "City", "State", "Type", "Location", "Date", "Day of Week", "Time Received", "Time Dispatched", "Time Arrived", "Time Cleared". I put up the data from January 2012 on Github in repo called SLO-Crime under my name, chadbunn.  I will be using CartoDB to geocode the data and then QGIS to perform the analysis.  I will probably make some maps in TileMill as well if I have enough time to. Any input would be appreciated because it is getting tedious! Thanks!

https://github.com/chadbunn/SLO-Crime/tree/master

Text sample below showing a single incident (plain text file, each section composed of a dispatch header summary followed by type / address / officer and other information):

===============================================================================
120101001 01/01/12 Received:00:05 Dispatched:00:06 Arrived:00:12 Cleared:00:53
===============================================================================
Type: MENTAL SUBJ                                             Location:LZ1
As Observed:
      Mental Health

Addr: 1257 RICH, San Luis Obispo, CA             Clearance Code:Report Filed in

Responsible Officer: Hyman, J
Units: 4235  ,4241
 Des: incid#=120101001 AP/KAVANAGH,KOLLEEN 060455 166 PC clr:RTF oc:MENH call=1l


texts: How do I add a georeferenced JPEG with a JGW Worldfile to OpenLayers v5.3.0?

texts:
I have a JPEG file with an associated .JGW worldfile containing the following data:

0.0144218
0
0
-0.0144218
-65758.8265181
-130841.5565499

My Jpeg is 21177*14242 pixels resolution.

I'm currently working with the following code:

let extent = [0, 0, 0, 0,]; //I need the correct values!

let exampleProjection = new ol.proj.Projection({
code: 'orto-image',
units: 'pixels',
extent: extent //I need the correct values!
});

let exampleImageSource = new ol.source.ImageStatic({
url: "https://mywebsite.com/uploads/example.jpg",
projection: exampleProjection,
imageExtent: extent //I need the correct values!
});

let exampleImageLayer = new ol.layer.Image({
source: exampleImageSource
});

I would like to present this image over my OSM map. How can this be achieved?


texts: LineString in OpenLayers: Is it possible to draw segments in different colors?

texts:
I'm a newbie in OpenLayers. Can you help me with this:

Is it possible to assign color not to the whole linestring, but to different segments?


texts: Invert Deleting multiple fields with delete tool in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.1 SP1.

I have a feature (polygon) with nearly 60 fields in it. Now I want to delete all fields except 20(or a list of) and they are not one after one, no they are somewhere in the droplist. But to select for a tons around 40 feature shape files would be a huge searching and clicking time...because the update for every new loaded feature his droplist and i can start searching again and again.

Normally it must be go the opposite way...or?

I am looking for a tool or script by "inverse droping". From the web I did not get a useful solution.


texts: Are easements normally marked on county property maps?

texts: There is a main gas line easement on the edge of our land, but the gas company is concerned about where the garden is. Would the line be marked well enough on a county map to say where the easement is?


texts: Populate print template title text from a dojo textbox

texts:
I have a print widget that starts up great. It's based on the typical ArcGIS print widget sample.

However, I want to tweak this a bit. There is a "titleText" property in the layout options that I want to set from a text box that the user types a title into.

So, based on the dojo textbox value, I assign that to the titleText property.

A lot of the code I find is in AMD so it's really tough to decipher (as I am still starting out in JSAPI).

Thanks in advance for any help!


texts: Splitting concave polygons into convex polygons in QGIS

texts:
I want to optimize and decompose several concave polygon into convex ones but not too small or to complicated ones.Ideally cutting in places that would minimize the creation of convex polygons and in a way they are more close to being regular.
(I have selected some of theme just to be more clear)

I have tried (and they don't work):

Triangulation - Creating a lot of small polygons with Delaunay or other methods creates a lot of small polygons and it doesn't help me in my final objective.

Convex-hull - It changes the shape of my polygons and I don't want to change that because they contain information.

I need this because each of these polygons have information that I'm going to keep in points in a mesh to make thermic calculations after. I want to create this elementary polygons to then take the centroid of each one of them and create this mesh. So because of that it would be bizarre taking a centroid out of the polygon or taking the pole of inaccessibility and having a polygon that its almost "twisted".
An output for this selected polygons would be to cut them in concave pieces and then when we take their centroids would be something like these:


texts: Receiving unknown exception when trying to perform a coordinate transform with PyQGIS

texts:
I'm trying to create a script that will use the decimal degree coordinates of point features. To do this, I want to perform a coordinate transformation on features whose coordinates are in meters. In the code below, I get the active layer and store its coordinate system (UTM zone 18N), then choose an NAD83 geographic coordinate system (EPSG:5498), and construct a QgsCoordinateTransform. I then attempt various methods of executing the transformation, with bost QgsCoordinateTransform.transform(QgsPointXY) and QgsGeometry.transform(QgsCoordinateTransform), but all of them result in Exception: unknown.

Python Console

li = iface.activeLayer()

liearth = ifacel.activeLayer().crs()

flatearth = QgsCoordinateReferenceSystem(5498)

lipts = [f for f in li.getFeatures()]

toflat = QgsCoordinateTransform(liearth, flatearth, QgsProject.instance())

for f in lipts:
  gPoint = toflat.transform(f.geometry().asPoint())
  print(gPoint)
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 2, in <module>
Exception: unknown

fPts = [f.geometry().asPoint() for f in lipts]

type(fPts[0])
<class 'qgis._core.QgsPointXY'>

QgsCoordinateTransform(li.crs(), flatearth, QgsProject.instance()).transform(fPts[0])
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
Exception: unknown
type(fPts[0])
<class 'qgis._core.QgsPointXY'>

fGeo = [f.geometry() for f in lipts]
type(fGeo)
<class 'list'>
type(fGeo[0])
<class 'qgis._core.QgsGeometry'>

newpoint = fGeo[0]

newpoint
<qgis._core.QgsGeometry object at 0x114d5edc8>

newpoint.asPoint()
(592932,4.52942e+6)

newpoint.transform(QgsCoordinateTransform(li.crs(), flatearth, QgsProject.instance()))
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
Exception: unknown

The only helpful detail that I've been able to find online on this is someone who also received an unknown exception a few months ago a little before QGIS 3 was released. However, hey were able to fix the issue using a QgsRectangle, which I cannot do.


texts:
Currently I am working with some satellite images (grib files) in pyqgis. As I want to work only with the cells in my study area, I need to clip the raster by a shapefile. I know I can do this by Raster tools in QGIS, but since I have a lot of rasters to work with, I wanted to use python and save time.

Perhaps somebody knows the pyqgis commands to do this?

I am new to QGIS and python.


texts: Spatial reference data from all layers in MXD

texts:
I have a MXD with many group which has many layers. I have tried to get coordinate system data from all layers with the code from Iterating through layers in mxd and describing coordinate system. So far no luck.

I have tried with this code. I also tried with comments on the given answer but no luck. This is my code:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]
layerlist = arcpy.mapping.ListLayers
for layer in layerlist(mxd, '', df): # Loop through layers
    sr = arcpy.Describe(layer).spatialReference
    if spatialReference.type == "Unknown":
        print layer + " has an unknown spatial reference\n"
    else:
        print layer.name + sr.name + "\n"
del mxd

And this is error I get.

Runtime error
Traceback (most recent call last):
  "<string>", line 6, in <module>
 AttributeError: DescribeData: Method spatialReference does not exist

And if I try this code, from official ESRI site, nothing happens:

import arcpy
arcpy.env.workspace = "C:/PATH_TO_FOLDER/"
feature_classes = arcpy.ListFeatureClasses()
for fc in feature_classes:
    spatial_ref = arcpy.Describe(fc).spatialReference
    if spatial_ref.name == "Unknown":
        print("{0} has an unknown spatial reference".format(fc))
    else:
        print("{0} : {1}".format(fc, spatial_ref.name))
del mxd


texts: Replacing Polygon layer with extent of another polygon in ArcGIS Desktop?

texts:
I need to insert a multi-polygon feature class into an existing vector layer (replacing the entire area within the former's footprint).

Is there a better way than just DELETE then MERGE polygons (in ArcMap?)


texts: How to parse EXIF GPS information to lat,lng decimal numbers

texts:
"EXIF_GPSAltitude": "(220.279)",
"EXIF_GPSLatitude": "(55) (40.711) (0)",
"EXIF_GPSLatitudeRef": "N",
"EXIF_GPSLongitude": "(8) (30.2282) (0)",

How should i understand above example as degree, min, sec?

Some exif data that i extracted has metadata listed as above. Is these formats specified anywhere ? or should I expect alot of different formats of gps positions in EXIF?

I need to write a method that gives the lat,lng as decimal numbers based on exif data, and I am unsure about how many different formats I should expect to be able to parse(I will learn this over time) and this question is mostly about above examples. I know the degree is 55 and 8, but not sure if its 30 mins and 0.22*60 secs, and in that case, why is the 0 there?


texts:
When using Landsat 8 imagery to calculate NDVI, will the NIR and Red bands be affected by cloud cover?

In other words, will clouds produce a wrong NDVI calculation?