texts: QGIS doesn't recognize BETA2007.gsb for Transformation from GK3 to UTM

texts:
I want to transform data from GK3 (EPSG:31467) to UTM(EPSG:25832) unsing the transformation grid BETA2007.gsb. BETA2007.gsb is already in the  C:\Program Files\QGIS 3.10\share\proj directory.

I should be able to select it as my default transformation method under Settings->Options->CRS->Default Datum Transformations. However, It doesn't appear as an option. It should be in the list shouldn't it?

Does anyone know how to add BETA2007.gsb into the list of default datum transformation methods? I tried it using version 3.10; 3.16 and 3.18. 1 year ago it worked under QGIS 4.X. It is strange that it doesn't work anymore or I'm doing sth. wrong now.


texts: Working with a Trimble GPS Pathfinder Pro XH and Juno SB to collect high quality GIS point data

texts:
Dug an unused, uncharged Trimble Pathfinder Pro XH (XH) out of my companies junk closet. I have learned that this device is supposed to be able to collection GPS points with 30 cm resolution, exceeding the accuracy of the Bad Elf GPS Pro+ devices that we are using in the field currently (paired to an iPad through ESRI Collector App).

I have started trying to figure out how to use the XH, but I have run into several road blocks:

The XH cannot pair with an iPad through the collector App / iOS
The Juno SB and XH seem to connect over bluetooth, but there seems to be no easy way to confirm this - and as an answer below points out, this connection is spotty and weak.
Even if I get them to communicate, through bluetooth or serial cable or something else, and register a signal to collect a point (which hasn't happened successfully yet), I still have to figure out how to move the point from either the built in GPS software (GPS Controller) or an ancient ArcPad installation (8.0) into ArcMap Desktop 10.6

So, the question I am focusing on right now is:

For those who are using a Trimble GPS Pathfinder Pro XH  and Juno SB
currently to collect GPS points, how do you bring data from this into
ArcMap or another GIS system?


texts: Google-Earth-Engine + Jupyter : filterBounds not working in Landsat imageCollection

texts:
I have a simple problem of not being able to select a single Landsat image using filterBounds within Jupyter environment. It used to work just fine in Earth Engine's online interface. I don't know why it should fail now.

import ee
ee.Initialize()

lat = 45
lon = -100
Radius = 1000

roi = ee.Geometry.Point([lat, lon]).buffer(Radius)
point = ee.Geometry.Point (lat, lon)
start = ee.Date('2017-01-01')
finish = ee.Date('2018-12-01')
L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1")

Time_Set = ee.ImageCollection(
  L8.filterDate(start, finish)
    .filter(ee.Filter.lt('CLOUD_COVER', 2))
)
print('Number of images in Time Set: ', Time_Set.size().getInfo())

Space_Set = ee.ImageCollection(
    L8.filterBounds(ee.Geometry.Point(lat, lon))
)

print('Number of images in Space Set: ', Space_Set.size().getInfo())

Printed results:

Number of images in Time Set:  64343
Number of images in Space Set:  0


texts:
I am trying to set a layer's style using pyQGIS. My primary idea was to set the style from a string, but unfortunately I can only do it from a file. This code works as a charm:

        filename= QDir.tempPath() + "/style.sld"
        f = open(filename,'w')

        if os.path.exists(filename):

            f.write(style)
            f.close()
            ret=layer.loadSldStyle(filename)

            if ret:
                self.log(logging.DEBUG, "Style loaded")
            else:
                self.log(logging.DEBUG, "Style not loaded: " + errorMsg)

            os.remove(filename)

As you can see, I am writing the style string into a temporary file, loading the style and then removing it. What I would like to do instead, is to pass the string straight to the readSld function, which is wrapped in loadSldStyle:

        d = QDomDocument()
        ret=d.setContent(style)

        if ret:
            self.log(logging.DEBUG, "Document created")
        else:
            self.log(logging.DEBUG, "Document not created")

        self.log(logging.DEBUG, "content: " + d.toString())

        myRoot=d.firstChildElement
         ("StyledLayerDescriptor")

        if myRoot is None:
            self.log(logging.ERROR,  u"Could not get "
                "root element")
            return

        namedLayerElem=myRoot.firstChildElement
         ("namedLayerElem")

        if namedLayerElem is None:
            self.log(logging.ERROR,  u"Could not find "
                "NamedLayer")
            return

        errorMsg = str("")
        ret=layer.readSld(namedLayerElem,errorMsg)

Please note that this is very similar to the C++ code for the method loadSldStyle on the QgsMapLayer class implementation

The QDomDocument contains the correct style, and the myRoot and namedLayerElem nodes are correctly created. Nevertheless, readSld returns false, and worst of all, the error string is empty so I have no idea what is going on!
Any advice on what might be the problem? I know that this function has been ported to sip.


texts: Extract Values To Points using ArcPy with ArcGIS Pro gives arcgisscripting.ExecuteError: ERROR 010555: Internal error in Sample

texts:
I have built an ArcPy tool that works using Python 2.7. I am trying to use the same tool but now using Python 3. When trying to perform ExtractValuesToPoints of ArcPy I keep getting the following error:

  arcgisscripting.ExecuteError: ERROR 010555: Internal error in Sample:
  %s. ERROR 002834: The geometry operation was cancelled. Workspace or
  data source is read only. Workspace or data source is read only. The
  table was not found. [VAT_mixed_ht_raster_with_blocks] The table was
  not found. [VAT_mixed_ht_raster_with_blocks] ERROR 002834: The
  geometry operation was cancelled. Failed to execute
  (ExtractValuesToPoints).

Executing with Python 2.7, the tool runs fine.

code:

arcpy.gp.ExtractValuesToPoints(pt_file, ht_raster, observers)

Does it have anything to do with gp?


texts: Find and Replace shapefile source script not working

texts:
I'm very new to python so please forgive me if this is a simple error.
I have quite a few MXD's in which each have a few layers referencing the same source. However I have since moved this data to a new folder.
I have written a script to automatically refresh the similar layers.

The script runs through, however when I check to see if its worked, none of the layers are refreshed. I'm not sure whats wrong.

import arcpy, os

#workspace to search for MXDs
Workspace = r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\MXD"
arcpy.env.workspace = Workspace

#list map documents in folder
mxdList = arcpy.ListFiles("*.mxd")
print mxdList

#set new link for specific shapefile in each MXD in list.
for file in mxdList:
#set map document to change
    filePath = os.path.join(Workspace, file)
    mxd = arcpy.mapping.MapDocument(filePath)

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        if lyr.dataSource == r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\WP1\MMO_South_Plan_Areas.shp":
            lyr.findAndReplaceWorkspacePath(r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\WP1", r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\SHP")
    if lyr.dataSource == r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\WP1\Land.shp":
            lyr.findAndReplaceWorkspacePath(r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\WP1", r"M:\TEAM_GIS\EDKM\Projects\mmo1065_south_habitat_mapping\20140708_Update_3_QA\3_Working_GIS_Files\SHP")

arcpy.RefreshTOC()
arcpy.RefreshActiveView()
mxd.save()
print "Success"


texts: I there a way to calculate the accuracy of a classified image (users-, producers accuracy, Kappa) in R, assuming I have a two columns matrix. On of the columns would contain the values that came out with the classification, the other one my observed/checked values?


texts: Checking given a coordinate (lat,long) if is in a raster image file (tif file) and how to extract it with NXN window

texts:
I am researching to complete a project.
The aim is by given a coordinate I can extract the image among the raster files.
How do I check if the coordinate that I have for example (51.3334198, 3.2973934) is in a raster image file - k_01.tif? If this coordinate is indeed in k_01.tif, how do I extract a small part of if i.e NxN window.
My code:
import rasterio src = rasterio.open('k_01.tif')

src.bound

>BoundingBox(left=145000.0, bottom=238000.0, right=162000.0, top=247000.0)

src.crs

>CRS.from_epsg(31370)

Thanks @snowman2 for giving the link Extracting data from a raster.
I followed the advise and my code as follow:
note: I used another .tif and another longitude and latitude

import rasterio
from pyproj import Transformer

lon = 4.2811439
lat = 51.2010511

path = './k15.tif'

with rasterio.open(path) as rds:
    # convert coordinate to raster projection
    transformer = Transformer.from_crs("EPSG:4326", rds.crs, always_xy=True)
    xx, yy = transformer.transform(lon, lat)

    # get value from grid
    value = list(rds.sample([(xx, yy)]))[0]

value
> array([11.61], dtype=float32)

xx , yy

>(143876.25575738482, 210123.58613220416)

src.bounds

> BoundingBox(left=130000.0, bottom=198000.0, right=162000.0, top=218000.0)

From the result above I understand that xx,yy  is the conversion of the latitude and longitude and it is indeed in the boundary of the tif. What is the value 11.6 for?
My second part of the question is how to extract the latitude and longitude with NXN window ? As shown as picture below.

And crop it to just as follow?

thanks in advance.


texts: ArcPy Label Class object returning unexpected value for expression property?

texts:
import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd, "TestFile_LineToPoints13_Prod"):
    if lyr.supports("LABELCLASSES"):
        for lblclass in lyr.labelClasses:
            lblclass.expression = " [OBJECTID]"
            if lblClass.showClassLabels:
                    print "    Class Name:  " + lblClass.className
                    print "    Expression:  " + lblClass.expression
mxd.save()
del mxd

Prints:

  Class Name:  Default

  Expression:  [Shape_Length]

The label for that layer is even set at ObjectID and is currently showing the labels... Why does it keep responding Shape_Length?


texts: PyQgis: Extract SVG-Symbols and SVG-Paths of all Layers of my QGIS Project

texts:
I want to extract the SVG-Paths and SVG-Symbols for all Layers in my QGIS Project.

The purpose is to open the SVG files via a plugin (so via PyQgis) and I don't know how to do so.
I checked the QGIS Pyton API but did not find a solution (have to admit that I'm quite new to PyQgis).

Would be glad if someone could give me a hint.


texts:
I am creating ArcGIS Online map and having a functionality to search Zip Code using the USA Zip Code 2014 pre-existing layer.
I further need to extend the functionality so that user can provide Zip Code & radius and can see visual representation of specified radius (from center of Zip Code)

Has anybody come across a similar scenario?


texts: SQL: How to union polygon geometries from a single table that have the same id and certain type?

texts:
Essentially, the table looks like this (geometries are mostly unique, and shouldn't overlap, but may overlap a little):

id         sub-id    type         Geometry

123456     1         ordinary     POLYGON(..)

123456     2         ordinary     POLYGON(..)

123456     3         shoulder     POLYGON(..)

456789     1         shoulder     POLYGON(..)

456789     2         ordinary     POLYGON(..)

456789     3         ordinary     POLYGON(..)

456789     4         ordinary     POLYGON(..)

The ST_union should combine polygons from rows 1 and 2 (they share same id and have certain desired type: ordinary. I don't care about other types).

Then, it should keep going and the same ST_Union should combine polygons from the last three rows shown (again, they share same id and have certain desired type: ordinary. I don't care about other types).

The code would then do an analysis using those aggregate polygons like, say, finding if there is an intersection.

There are perhaps hundreds of thousands of sub-ids so I'm hoping that merging the geometries of the sub-ids (by id and the 'ordinary' type) and then running the intersection analysis later in the code will make increase the efficiency of the query.

I don't want to/shouldn't change the original database, and I can't create a new database or table to be saved on the server.


texts: Output raster from geostatistical layer to match template raster: why is there an extra row of NoData cells?

texts:
I have a geostatistical layer created using ordinary kriging and the Create Geostatistical Layer tool in Model Builder. I am trying to output this layer to a raster that matches a target raster using the GA Layer to Grid tool.  I have used the target raster to set the Snap Raster, Cell Size, and Processing Extent environment variables and various logical combinations of these variables.  However, the raster that is output has one additional row at the top when compared to the target raster.  This row is filled with NoData values.

Differences are highlighted in bold italics

The target raster has these dimensions:

  Columns_and_Rows: 110, 166
  Cell_Size: 0.11559105, 0.11559105
  Extent:

    Top     49.359975678
    Left  -129.652945807
    Right -116.937930333
    Bottom  30.1718614184

  Spatial_Reference: CS_WGS_1984

The output raster has these:

  Columns_and_Rows: 110, 167
  Cell_Size: 0.11559105, 0.11559105
  Extent:

    Top     49.4755667278
     Left  -129.652945807
     Right -116.937930333
     Bottom  30.1718614184

  Spatial_Reference: CS_WGS_1984

Why does my raster have an extra row (and therefore the wrong Extent) and how can I fix this?


texts: OpenLayers and Google Maps when panning tiles not in sync

texts:
I want to overlay WMS tiles on Google maps using the following scales
1:287 / 1:573 / 1:1146 / 1:2292  (layer preview in geoserver)

I set a bounding box and a centroid for a specific area to set the google map as base layer.

However as I zoom or move the mouse my overlayed tile moves along, losing focus (when panning tiles not in sync).

Is there a way to resolve this?

(i am using openlayers vs 2.1/geoserver 2.1.4)


texts: Is there a way to figure out what files an MXD points to?

texts:
Let's say I load up two layers, and then save an .mxd file.  Now I want to get the path names of my layers from the .mxd file.  Is there a way to do this in ArcEngine?

I'm using ArcEngine 10 with VS2010 & C#


texts:
I am unable to remove the legend from a CartoDB layer. I have passed legends: false in createLayer options but with no luck. Below the code i used :

cartodb.createLayer(map, BikeTrails)

                .on('done', function(layer) {legends:false

                layerBikeTrails = layer;

L.control.layers(null, {
    'Dedicated Bike Lanes': layerBikeLanes,
    'Bike Share Toronto': layerBikeShare,
    'Bike Stores': layerBikeStores,
    'Access to Multipurpose Trails': layerBikeTrails
})
.addTo(map);


texts: Update fields automatically after adding new geometry

texts:
I have two layers, one for "States" (Polygon) and the other for "Important_Cities_location" (Point). "States" feature class is completed with all attributes. I want to update the state name field automatically in "Important_Cities_Location" when I place them in the respective "States" spatially.

Can anyone help?


texts: Creating a python string to be passed as anSQL Expresion(forqgsExpression) from a dict

texts:
Let's say I have the following dictionary:

  foodict = {'a': ['a1','a2'], 'b':['b1','b2'], c:['c1']}

I want to generate the following python string that will be passed as an SQL expression:

'"a" = 'a1' AND "b" = 'b1' AND "c" = 'c1' OR "a" = 'a1' AND "b" = 'b2' AND "c" = 'c1' OR "a" = 'a2' AND "b" = 'b1' AND "c" = 'c1' OR "a" = 'a2' AND "b" = 'b2' AND "c" = 'c1''

Except the opening and closing quotes, the single/double quotes matters since it is intended as an sql expression being passed to pyqgis.

Example :

def extract_by_field_category(inputData, selectionDict, outputData):
lyr = QgsVectorLayer(shp_path,'lyr','ogr')
selection = lyr.getFeatures(QgsFeatureRequest().setFilterExpression(SQLEXPRESSION))

The dict is an hardcoded parameter in a function call, any suggestion on how to improve it is welcome to. The idea of a dict as argument /parameter seems convenient to reduce the typing needed in the function call.


texts: PostGIS intersection betwen polylines and buffer around points (polygons)

texts:
I have a PostGIS table of polylines and I need to intersect them with a table of polygons (conveniently buffers around points). In real  world polylines are cables, and points are poles. I need to know for each cable what poles are along this cable in order
What is very important is the resulting intersections have to be in order as I read the polyline from one end to the other.


texts: Why does Field Calculator put in all these 9's?

texts:
I’m trying to update some GPS point coordinates in an attribute table using QGIS Field Calculator. However, regardless of whether I do it by creating a new field or updating an existing one, I get about 13 decimal places, the last 5 of which are “9” for every feature. (The original points only had 6 decimal places.) I doubt that this is really a coincidence! It happens even if I stipulate a lower Output field width and Precision. My expression is $y, for example, for latitude.

Any idea what’s going on here?


texts: Reprojecting over the Pacific Ocean with EPSG:3832 gives artefacts

texts:
I am trying to center the Natuarl Earth 10m on the Pacific. I have followed instructions from QGIS Tutorial: How to Change the Map Projection to be Centered Over the Pacific Ocean.
When I change projection to EPSG:3832 the map in QGIS becomes fragmented:

How can I fix this? Do I need to define my own CRS entry?
I have seen other replies here Displaying world country shapefiles centered on Pacific Ocean using Robinson or Miller Cylindrical projection in QGIS? but that is 7 years old.

This is an edit: Interestingly open street map seems to reproject to 3832 ok. Slow and blurry though.


texts: Why do Landsat images sometimes form triangles instead of 4 sided polygons to cover areas?

texts:
I put several Landsat8 images into ArcMap 10.7. I don't get the same kind of coverage (square area) you see in the first picture from the Earth Explorer site.  I instead get picture 2: corners, missing, triangular shapes for some etc.
Picture 1:

Picture 2:


texts: ArcPy 10.0 iterate through shapefile1 clip shapefile2 to shapefile1's polygons

texts:
I know that this question has been answered several times before, but for some reason I can't get this script to run... So I was hoping someone could point out where I'm going wrong.

I have two shapefiles, shapefile 1 has polygon boundaries and shapefile 2 has a single polygon that overlaps the firsts polygonal boundaries.

I would like to clip shapefile 2 according to each polygon of shapefile 1.

I have seen several sites and examples, but for some reason I must be overlooking a detail in getting this to work

Some of the examples I've seen are:
How to Batch Clip in ArcGIS Desktop 10 using Python/ArcPy?

Batch clip feature class to preserve area field

Clipping buffers to census tracts

import arcpy
arcpy.env.workspace = "path\\to\\workspace\\MaintenanceDivisions"
arcpy.env.overwriteOutput=True
#variables
md = "MaintenanceDivisions.shp"
out = arcpy.CreateUniqueName( "clipped.shp" )
mfp = "dissolved.shp"
try:
    #clip polygons from selected area
    rows = arcpy.SearchCursor(md)
    row = rows.next()
    for row in rows:
        feat = row.Shape
        arcpy.Clip_analysis([mfpd], feat, out, "")
        row = rows.next()
except arcpy.ExecuteError as e:
    print(e)
del rows


texts: How to merge SeaWiFS and MODIS data or use other merged products?

texts:
I hope to produce a time-series for sea surface temperature (SST) and chlorophyll-a (Chl-a) for use in MaxEnt without the use of arcgis, preferring Python tools to others. I am new to satellite data processing, and there appears to be quite a lot to consider.

What must be taken into account for merging SeaWiFS and MODIS datasets (e.g. discrepancies between the data, projections, etc.)?

Are there other products available that provide such merged data that are generally recommended?

I have found the following:

NOAA - CoRTAD (SST)
ESA - MEDSPIRATION (SST)
ESA - GlobColour (SST & Chl-a)

Ultimately I need all of the data gridded to the same grid size/area, so I would to write a script where I could acquire the data, re-project if necessary, re-sample to the grid if necessary, and then save to desired format (ESRI's ASCII .asc) for running in the model.


texts: Some edges of raster switch to a value of 0 when clipped and pixels are shifted on QGIS

texts: I am trying to clip trees out of a raster layer using polygon I've traced around them. I've created one polygon that I've copied for each tree so that each clipped raster contains the same width and height in pixels. But for some reason, some pixels from the edge become 0 values after the clip, and the pixel will have shifted. I'm guessing it is because my mask layer doesn't cover some pixels completely around the edges. There is no way to snap a polygon to edges of raster pixels that I know of. So how could I avoid this issue.


texts: hillshade and slopeshade geotiffs fail to convert to new projection using gdalwarp

texts:
I'm trying to convert my three DEM products (colour, hillshade and slopeshade; the latter two generated using gdaldem) from EPSG:4326 -> EPSG:3857 using gdalwarp, via the command:

gdalwarp -of GTiff -s_srs epsg:4326 -t_srs epsg:3857 -co "tfw=yes" {src} {dst}

This works fine for the colour geotiff, but the resulting hillshade and slopeshade geotiffs now appear blank.

Sample output from gdalinfo for the slopeshade geotiff is below, with links to the original and converted files at the bottom.

gdalinfo before:

Driver: GTiff/GeoTIFF
Files: slopeshade_dem.tif
Size is 4557, 3577
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-121.460761599999998,38.583444700000001)
Pixel Size = (0.000000654597323,-0.000000654598826)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-121.4607616,  38.5834447) (121d27'38.74"W, 38d35' 0.40"N)
Lower Left  (-121.4607616,  38.5811032) (121d27'38.74"W, 38d34'51.97"N)
Upper Right (-121.4577786,  38.5834447) (121d27'28.00"W, 38d35' 0.40"N)
Lower Right (-121.4577786,  38.5811032) (121d27'28.00"W, 38d34'51.97"N)
Center      (-121.4592701,  38.5822740) (121d27'33.37"W, 38d34'56.19"N)
Band 1 Block=4557x1 Type=Byte, ColorInterp=Gray
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=4557x1 Type=Byte, ColorInterp=Alpha

gdalinfo after:

Driver: GTiff/GeoTIFF
Files: slopeshade_dem.tif
Size is 4088, 4105
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-13520950.132675198838115,4662178.034760829061270)
Pixel Size = (0.081230547006284,-0.081230547006284)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-13520950.133, 4662178.035) (121d27'38.74"W, 38d35' 0.40"N)
Lower Left  (-13520950.133, 4661844.583) (121d27'38.74"W, 38d34'51.97"N)
Upper Right (-13520618.062, 4662178.035) (121d27'28.00"W, 38d35' 0.40"N)
Lower Right (-13520618.062, 4661844.583) (121d27'28.00"W, 38d34'51.97"N)
Center      (-13520784.097, 4662011.309) (121d27'33.37"W, 38d34'56.19"N)
Band 1 Block=4088x1 Type=Byte, ColorInterp=Gray
Band 2 Block=4088x1 Type=Byte, ColorInterp=Undefined

Links to the slopeshade file before conversion (in EPSG:4326):

https://www.dropbox.com/s/9safqs1hkngnasf/slopeshade_dem.tif?dl=0

Links to the slopeshade file after conversion (in EPSG:3857):

https://www.dropbox.com/s/ebmgnbhrzog01ui/slopeshade_dem.tif?dl=0


texts: Finding if two polygons Intersect in Python?

texts:
I'm looking for an algorithm, a high level solution, or even a library which can help me determine if two polygons intersect, in Python.

I have the vertices of the two polygons (These are single part polygons without any holes) in two different arrays. The polygons are 2D (i.e. just X and Y coordinates)

I'll like to make an function which will return a boolean indicating whether these two polygons intersect.

Please note that I cannot use arcpy, or any arcgis components in this.

Can you suggest an algorithm or library for doing this?


texts: Getting features in order they are selected using PyQGIS?

texts:
Is there a way of establishing the order that a user selected features using QGIS and Python?

By example: we have a network of line features representing a road network, and the user defines a route from point A to point B. Currently, using layer.selectedFeatures() returns the correct features, but in the order they were drawn rather than the order they were selected.

I've tried a few things to do this manually:

QgsMapToolEmitPoint(canvas) obviously doesn't work - I want the feature not the point clicked. But I guess in worst case I could hack a spatial query to do a look up of the closest feature. Not ideal.
QgsMapToolSelect(canvas) doesn't seem to be available as part of the Python accessible GUI or CORE, although I can find mention of it existing elsewhere.


texts:
I have tried to convert a raster to points using this syntax (r.to.vect):

r.to.vect input=g_kelowna_merged@kelowna_mapset output=g_kelowna_points feature=point

I have tried using the GRASS interface in QGIS and also directly in GRASS GUI 6.4.3 and 6.4.4.

What happens is that the points get created and show in the Map name box but there is nothing in the Layer box and QGIS says "No layers available in this map". Looks like I created a map, not a layer. Any suggestions for getting those points?


texts: Postgis List of adjacent LineStrings ordered sequentially from one end of the MultiLineString to the other

texts:
I have LineStrings that represent a road network/address centerline that are segmented at intersections, etc.  I have determined how to query collections of LineStrings that are adjacent to each other and have the same street name. Now I am hoping to be able to list the sub-segment IDs sequentially based on adjacency, regardless of line direction of the sub segments since they do not all “flow” in the same direction.  So that they are listed in order from one end to the other. How do I implement this?  I know postgis_topology is an entire suite of topology related functions, but I’m hoping there is a simple way to compare coordinates.

WITH
     endpoints AS (SELECT ST_Collect(ST_StartPoint(shape), ST_EndPoint(shape)) AS geom
                   FROM raod_segments),--convert segments down to start and end points
     clusters  AS (SELECT unnest(ST_ClusterWithin(geom, 1e-8)) AS geom
                   FROM endpoints), --Unsplit/merge adjacent roads
     clusters_with_ids AS (SELECT row_number() OVER () AS cid,
                                  ST_CollectionHomogenize(geom) AS geom
                           FROM clusters) --add parent IDs to merged segements
SELECT
    ST_Collect(raod_segments.shape) AS unsplit_geom,  ---create merged geometry as a collection of segments
    string_agg(raod_segments.objectid::text, ' | '
        --order by ST_YMin(ST_StartPoint(shape)) – HOW TO ORDER BY ADJACENCY STARTING FROM ONE END
        ) as child_segment_objectid,   -- list all child segments that make up the parent road
    fullname,
    COUNT(fullname) as split_count
FROM raod_segments
LEFT JOIN clusters_with_ids ON ST_Intersects(raod_segments.shape, clusters_with_ids.geom)
GROUP BY cid, fullname  --group the segments by parent id and unique street name.
;

Background:  These are road segments that are missing address ranges.  I am exploring methods to auto-populate theoretical address ranges for each segment in order to support geocoding by interpolation. What terms should I be googling? Pitfalls associated the proposed process.


texts:
I am a complete beginner in OSM and osmosis, and I am trying to build a stripped down version of a local map, in order to add nodes of public lights that would be used for maintenance. I have used the following command line in order to keep only streets and roads, but there are still plenty of one way signs (arrows)  cluttering the narrow streets in the center of the village, and I would like to get rid of them. I could not find any way to eliminate them (or to make them invisible in the output) I am also trying to get rid of the names of the streets, since most tend to consume too much space in some areas where I will need to add nodes.

osmosis --read-xml file="Source-Map.osm.xml" --tag-filter accept-ways highway=* --tag-filter reject-ways highway=track --tag-filter reject-nodes highway=bus_stop --tag-filter reject-nodes highway=one_way --tag-filter reject-nodes barrier=* --tag-filter reject-nodes highway=untagged --tag-filter reject-relations --used-node --write-xml file="FilteredMap.osm"


texts: Find all points near a point within a radius MySql 8.0 (SRID 4326)

texts:
This question is targeted to MySQL 8.0 with its new spatial features.

Given the following table:

CREATE TABLE `places` (
  `id` INT(11) UNSIGNED NOT NULL,
  `location` POINT NOT NULL SRID 4326,
  PRIMARY KEY (`id`),
  SPATIAL KEY `location_idx` (`location`)
);

I want to retrieve all the places near to a Point (let's call it MY_POINT) within a certain distance in meters (let's call it MY_DISTANCE), I also want to retrieve the distance within that point, and I want to do it in the most optimal way, e.g: triggering a spatial index with a WHERE clause.

To get the distance, I must use the ST_Distance function to calculate the distance between MY_POINT and each of the query results;
To test if the point is within MY_DISTANCE, I must create a circle around MY_POINT with MY_DISTANCE as a radius and then check if MY_POINT intersects the polygon.

This is the query I've made so far:

SELECT
  id,ST_Distance(MY_POINT, location) AS distance
FROM
  places
WHERE
  ST_Intersects(ST_SRID(ST_Buffer(MY_POINT, MY_DISTANCE / 111194.927), 4326), location) -- 111194.927 = degree of latitude in meters.
HAVING
  distance <= MY_DISTANCE -- This fixes the inaccurate ST_BUFFER polygon in some cases.
ORDER BY
  distance
ASC;

The problem is ST_BUFFER is making a weird and inaccurate circle (nothing to do with segments, rather I think it's because I calculate wrong the radius) and I'm missing some results.

Can anyone throw some light to this?


texts: Looking for intersection of polygons with a tolerance in QGIS

texts:
I have two QGIS layers I created from shape files.  I wanted to find which regions from one layer intersect with districts from the other layer, and the intersect tool does that perfectly.  I find however that in part of my map, the regions and districts coincide, but due to either rounding or tiny errors in my data, or something, there are extra intersections.  For example a region along the border of one district appears to (incorrectly) intrude a tiny bit into the next district.

I am looking for a way to create a list of intersections, but not count intersections where the area of the intersection is less than some tiny amount.  (I would be okay expressing the tolerance in any units, either as an area of as a distance)


texts: Annotating satellite images in the field

texts:
I have high res satellite imagery of my remote field location through Planet. I am working with indigenous communities on participant action mapping and would like to be able to annotate my images in the field. I have a microsoft surface pro and would like to be able to draw polygons on the images, as well as add information to associate with each polygon such as, in the example of a farm site, what crops are being grown there, who works there, what methods are being used.

I can run GIS on this machine but I don't know if I can draw polygons with a stylus onto the map in GIS or add this data using GIS. Does anyone know if this is possible and if so how to do it? Or do you have any ideas for a different software that I can use to accomplish this?


texts: QGIS Slope Analysis produces strange values

texts:
I am performing a slope analysis (Raster>Analysis>Slope) of a US Geological Survey DEM.
I am calculating the slope of a raster DEM saved as a TIF. I set up the slope analysis as per the picture below in EPSG 4326. I tick on calculate percentage and compute edges and leave all other variables as standard. When I calculate the slope in percentage I get really odd values being e to the 6 as shown in the layer panels in the second item down, rather than between 0 and 90%. In the map the layer looks correct though and this is a hilly area. How do I make it so the values appear between 0 and 90?
Similarly when slope analysis is performed in degrees (not ticking the percentage option it returns nearly all cells as 90 degrees which isn't correct (layer 2 in layers panel).
How do I solve these problems?


texts: Making Point Z feature class out of Point feature class?

texts:
I have a point shapefile where i need to transfer the records into a geodatabase feature class that is Point Z.  I have found a number of processes to go the other way, removing Z value and make it a simple point, but I want to go the other way.  I know I do not have Z values so wondering if this is possible and the points would just get a Z value of 0, or if Arcmap will just continue to give me errors when I try to convert.  I tried enabling Z values in the geoprocessing environment menu and then did a Feature Class to Feature Class geoprocess, but it created the new feature class without Z values still.

Any suggestions?


texts: GeoServer fails to read a JPEG compressed TIFF with alpha band

texts:
I have been trying to publish a GeoTIFF (with JPEG compression + alpha band) but couldn't succeed. The TIFF is generated using this GDAL command :

gdalwarp -t_srs EPSG:4326 -r bilinear -dstalpha -of GTiff -co COMPRESS=JPEG -co JPEG_QUALITY=25 -co TILED=YES input.vrt output.tif

I am using the dstalpha option of gdalwarp in order to avoid some black seamlines on the borders of images.

Unfortunately, it seems that GeoServer doesn't handle Alpha band with JPEG compression, I get this kind of exception :

  Caused by: javax.imageio.IIOException: Unsupported Image Type
  geoserver
  |       at java.desktop/com.sun.imageio.plugins.jpeg.JPEGImageReader.readInternal(JPEGImageReader.java:1139)
  |       at java.desktop/com.sun.imageio.plugins.jpeg.JPEGImageReader.read(JPEGImageReader.java:1110)
  geoserver
  |       at it.geosolutions.imageioimpl.plugins.tiff.TIFFJPEGDecompressor.decodeRaw(TIFFJPEGDecompressor.java:282)
  |       at it.geosolutions.imageio.plugins.tiff.TIFFDecompressor.decode(TIFFDecompressor.java:2637)
  |       at it.geosolutions.imageioimpl.plugins.tiff.TIFFImageReader.decodeTile(TIFFImageReader.java:1764)
  |       at it.geosolutions.imageioimpl.plugins.tiff.TIFFImageReader.read(TIFFImageReader.java:2065)
  |       at com.sun.media.jai.imageioimpl.ImageReadOpImage.computeTile(ImageReadOpImage.java:697)

When doing gdalinfo on my file, everything seems fine :

Driver: GTiff/GeoTIFF
Files: subset_nz.tif
Size is 4404, 4404
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (175.499992707622397,-37.099986337232799)
Pixel Size = (0.000022707000000,-0.000022707000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 175.4999927, -37.0999863) (175d29'59.97"E, 37d 5'59.95"S)
Lower Left  ( 175.4999927, -37.1999880) (175d29'59.97"E, 37d11'59.96"S)
Upper Right ( 175.5999943, -37.0999863) (175d35'59.98"E, 37d 5'59.95"S)
Lower Right ( 175.5999943, -37.1999880) (175d35'59.98"E, 37d11'59.96"S)
Center      ( 175.5499935, -37.1499872) (175d32'59.98"E, 37d 8'59.95"S)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha

I am running GeoServer 2.16.2, is this a known issue or is there anything wrong with my files/config?


texts: Understanding private "with link" visualizations in Carto?

texts:
We are creating web dashboards for a small number of customers (~5) in which each customer can view private business-related statistics.

We want to embed some visualizations inside the dashboards. The visualizations are static, i.e. the customers will not be able to change the "SQL query" (by "static" we do not mean "static image" - zooming, dragging and mouse-hovering should still work).

The data and the visualizations should be kept private as they may have commercially sensitive value.

Password protected visualizations will be annoying to the costumers - they already entered a password when they logged in to the dashboard.

Is the "With link" privacy mode what we need? i.e. simply create an iframe to embed the visualization, and assume that guessing the url is practically impossible.

What happens if we accidentally switch to a "Public" mode?

Where exactly is the visualization published?


texts: Count number of features in a geopackage and write to a file

texts:
I'm tasked to create a simple catalogue of 2000+ geopackages, and hope there is a way to do that with QGIS.
The idea is to batch run over all these geopackages and write the following to a CSV file;
file name, number of features, geometry type, srid, bbox
for eg, a CSV entry like this;
countour.gpkg, 1876, multipolyline, 4326, bbox
..
..
Did anyone had such a requirement in the past? It will be great if this can be done using QGIS processing toolbox as a batch process.
EDIT
There are 2000+ gpkg with one table in each file. Add to that, each geopackage has strictly one table with one geometry type only..


texts: How to get the feature location in OpenLayers v3

texts:
After I am using the drawing interaction or modify interaction I have a listener that returns the feature that modified or added.

Example

draw.on('drawend', function (event) {
    // get the feature
    var feature = event.element;
    // ...listen for changes on it
    logStatus(feature.getId());
});

I know how to get the id, but I need the location (lon lat) of the feature because I need to save it to the DB, how can I do it? I didn't find it in the API.


texts: NAVD88 not in Vertical Reference drop down list of Drone2Map

texts:
I use the Propeller company's Aeropoint GPS tiles to collect ground control for drone flights.
When trying to import my ground control points, the survey report provides the coordinate reference system in:
Vertical coordinate reference system: NAVD88 height
Vertical EPSG code: 5703
Geoid model: geoid12b
However, Drone2Map has a very limited selection of Vertical Reference options to choose from:

EGM84
EGM96
EGM2008
ELLIPSOID
HEIGHT ABOVE ELLIPSOID

Can anyone help me choose the correct option?


texts: unexpected output from raster calculator function

texts:
Working in ArcMap, I am trying to convert the logit values in raster into their actual probabilities (logistic regression,) using the following function:

1/(1+Exp(-"sn_topo_poly_log"))

Ive used this function before and it should return a raster with values from 0 to 1. Instead I get a raster with all values=0.

the input raster "sn_topo_poly_log," looks correct and has values from close to zero to -871361


texts:
How can I batch create curved linear features i.e. using 3 points in ArcPY?

The real question is, whether we can create curved linear feature through arcpy or we have to use ArcObjects! Because, to me, it seems we can not do that in ArcPy.

to create linear (non-curved) features we can use:

arcpy.Polyline(Array([arcpy.Point(*coords) for coords in feature])))

This question is a duplicate of Generating true curve elliptical polygons in file geodatabase using ArcPy? New comers, note that curves are only supported on 10.3 and above. see here


texts: How can I create an Access (.mdb) file containing a personal geodatabase?


texts: Cartographic Representation: how to use different rule for different polygon?

texts:
I create a file geodatabase to hold a polygon of the land boundary. What I would like to set up the RULE to separate colors. I had to split the polygons into 2 separate polygons and change the ID to 2. How do I add that to the Resprentation Symbology properties ? I already have set up for the north part but need to change the color the south part.

Thanks in advance !

here is the picture I want to change the colors.


texts: Visualizing and downloading Landsat 5 NDVI Composite images

texts:
I'm working with GEE (no so experienced), trying to download a series of Landsat 5 NDVI composite images. The thing is the code I have isn't showing me the images that I want to downlad afterwards.
Here is my code:
    var L5_NDVI: ImageCollection "Landsat 5 TM Collection 1 Tier 1 8-Day NDVI Composite"
    var Zone: Polygon, 4 vertices

    var imageVisParam ={
      min: 0.0,
      max: 1.0,
      palette: [
        'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
        '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
        '012E01', '011D01', '011301'
      ],
    };

    var L5_NDVI_d1 = L5_NDVI.select('NDVI');

    var date1 = ee.Date('2004-01-01');
    var date2 = ee.Date('2004-01-31');

    var L5_NDVI_d1 = L5_NDVI
      .filterBounds(Zona)
      .filterDate(date1, date2)
      .filterMetadata("CLOUD_COVER", "less_than", 10);
    print(L5_NDVI_d1);

    var mask = function(image){
      var cut = image.clip(Zone);
      return cut;
    };

    var L5clipped = L5_NDVI_d1.map(mask);

    Map.addLayer(L5clipped, imageVisParam, "landsat");

      Export.image.toDrive({
      image: L5clipped,
      description: 'Landsat_Jan04',
      folder: 'GoogleEngine',
      scale: 30,
      region: Zone
    });


texts: Is there mathmatical function/equation available for ArcGIS Network Analyst's Closest Facility Method?

texts: I would like to use this function in other statistical programming, ideally SAS where td = (whatever the function is to calculate total distance from centroid of county to facility). The goal is to leverage other demographic weights in drive time and distance.


texts: Convert relative coordinates into absolute ones having only one absolute coordinate in QGIS

texts:
I have a layer of points with relative altimetric heights from a survey carried out with a Differential GPS. I know the absolute coordinates of only one point (EPSG:4326 WGS 84 ), the other points have coordinates in an arbitrary system.

Is it possible to transform all other relative coordinates into absolute ones based on a single point known with QGIS?


texts: LAS .PTC classes into ESRI LAS to Mulitpoint

texts:
We have some Lidar data that has been classified based on visable features (i.e. trees, streets, buildings, etc.) These LAS data are used inside inside Microstation w/ TerraScan connected to a PTC file.

PTC file resembles this;
0       Unclassified
*   Zero length line    9   0   0
1       Unknown
*   Zero length line    9   0   0
2       Bare Earth
*   Zero length line    2   50  0

Is there a way to use this PTC file when symbolizing the resulting ESRI point cloud (ArcScene) following the LAS to Multipoint conversion. All the points are showing, but need to find out how to use the FC in the symbology inside ArcGIS.


texts:
I am experiencing a very annoying phenomenon at the moment when trying to plot an A0 PDF Plot and would be interested to hear if anyone has had a similar problem and whether there is a solution or workaround.

I have exported a map consisting of vector layers with raster background in A0 Format to pdf with ArcGIS 10.3.1. The exported PDF dokument looks fine when opened in PDF Reader. When I plot it with our HP T1300 plotter using the HPGL Driver I can see a straight diagonal line crossing through the whole map. This line does not exist in the sde Featureclass!

I used to have similar problems with ArcGIS 9.3.1 and a simple geometry check and repair would usually solve this. As this the data comes out of an MSSQL SDE the check/repair geometry tools is not available for use. Apparently the geometries are already repaired by importing into the sde.

Is this a common problem and are there any tried and tested solutions?


texts: How can / should sf spatial objects be translated (i.e. moved within same projection)?

texts:
Big fan of the sf (and previously, sp) spatial packages in R.

I see spTransform and st_transform as methods to reproject a spatial object into a new space. However, I have a case where I have a complicated spatial polygon object centered on a point. I'd love to not have to (inefficiently) recreate the spatial polygons from scratch at a new point (doing this for many census units centroids), but "recenter" the entire object iterative over each point. Any suggestions on how best to do this with sf (or worst case, going backward, sp) objects?


texts: How to solve Projection Problem at client side

texts:
I want to discuss a scenario here to all.

I have developed a Web mapping application in OpenLayers and GeoServer. The customer give us shapefiles with native projection SWREFF99_1330 (EPSG:3008).  i published this shapefile on Geoserver after changing its native projection from EPSG:3008 to EPSG:4326 by using ogr2ogr tool. The layer displayed perfectly over Google map layer in my OpenLayers Web Page.
Again, I publish this shapefile on GeoServer with its native projection (this time, without changing projection using ogr2ogr tool) and then i go to my page and figure out that the layer was not visible over Google map layer.

Now, my question to all members:

Is there any way to change the native projection from EPSG:3008 to EPSG:4326 at client side (in my OpenLayers web page code) ? so that the layer should display correctly over google base layer as in case 1.


texts: All house numbers in given zip and street

texts:
I'm a very beginner with overpass. I need to get all possible house numbers in area with given zip and a part of street.

For example.
All houses for: 10178 alex

The result should be or similar to:

Alex-Wedding-Straße, 10178 Berlin
Tunnel Alexanderplatz, 10178 Berlin
Alexanderplatz, 10178 Berlin
Alexanderplatz 1, 10178 Berlin
Kleine Alexanderstraße 1, 10178 Berlin
Alexanderstraße 1, 10178 Berlin
Alexanderplatz 2, 10178 Berlin

I tried the next query, but it returned no results.

area[postal_code=10178][boundary=postal_code]->.a;
way(area.a)["highway"~"^(primary|secondary|tertiary|residential|living_street)$"]["name"~"^alex.*"]->.street;
(
  node["addr:housenumber"](area.street);
  way["addr:housenumber"](area.street);
  relation["addr:housenumber"](area.street);
);
(._;>;);out;

What am I doing wrong?


texts: How to stop ArcGIS Desktop automatically creating mapped drives?

texts:
I have a MXD which I am using to make a PMF. I'm using ArcMap 10.1. I want my layers only to have UNC paths (starting with \\server...) as otherwise some users will get broken links if we have inconsistent mapped drives. Unfortunately, sometimes ArcMap 10.1 switches some layers to my mapped drive equivalent of the file path. I don't know why. When I try to change them all back to UNC paths in a wholesale fashion by disconnecting the relevant mapped drive so I can then switch one layer to the UNC path and let that propogate through all, ArcMap 10.1 decides instead to re establish the relevant mapped drive when starting up.

Anyone know a fix?


texts: What FOSS GIS has the option to link attribute tables to rasters

texts: In ArcGIS rasters can have attribute tables. What FOSS GIS offers a similar option?


texts: How can I change number format and alias in javascript API legend?

texts:
I want to change the alias and the length of a field that already exists in my GDB from the Javascript API. Is this possible?

I see that alias is a property of the Field class. Can I do something like featureLayer.field["fieldName"].alias = "fieldAlias"?

EDIT: I've changed the alias of the field in ArcMap and that persists in the feature service. But, I'm still wondering: How do I change the number format used in the dijit Legend? (I would like to display rounded numbers and percentages.)


texts: ogr: what is the advantage of using `geom.Clone()` over `geom`?

texts:
What is the advantage of using geom.Clone() over geom?

Simple example:

geom = link_feat.GetGeometryRef()
newfeature.SetGeometry(geom.Clone()) # or
newfeature.SetGeometry(geom)
layer.CreateFeature(newfeature)


texts: How to programmatically add an algorithm from external plugin to QGIS processing toolbox?

texts:
I'm writing a python plugin for QGIS and I'd like to know if it's possible to provide an interface to it from the processing toolbox.

I can see how the existing algorithms are loaded into processing via providers, but I don't see any obvious way of loading an algorithm without modifying the module. Perhaps it isn't possible and I need to submit a feature request, or attempt at coding it myself.

I was expecting something like:

import processing
from myplugin import mysweetalgorithm
processing.load_algorithm(mysweetalgorithm)


texts: How to get GeoJSON features using OpenLayers?

texts:
Can anybody explain how to get geojson features .Can somebody show way to get attributes from geojson text.

OpenLayers.ProxyHost = "/gt/proxy.cgi?url=";
     var lon = 5;
    var lat = 44;
    var zoom = 5;
    var map, layer;

    //function init(){
       var bounds = new OpenLayers.Bounds(
                -125051.196, 3967670.8507,
                1549596.1896, 4677667.0145
            );
            var options = {
                controls: [],
                maxExtent: bounds,
                maxResolution: 6541.59135,
                projection: "EPSG:23036",
                units: 'm'
            };
    map = new OpenLayers.Map('map', options);
    map.addControl(new OpenLayers.Control.PanZoomBar({
                position: new OpenLayers.Pixel(2,10)
            }));

            map.addControl(new OpenLayers.Control.Navigation());
            //map.addControl(new OpenLayers.Control.Scale($('scale')));
            map.addControl(new OpenLayers.Control.MousePosition({element: $('location')}));
             var geojson_format = new OpenLayers.Format.GeoJSON();
        var vl2 = new OpenLayers.Layer.Vector("",
        {isBaseLayer: true,
         styleMap: new OpenLayers.StyleMap({'default':{
                strokeColor: "#00FF00",
                strokeOpacity: 1,
                strokeWidth: 1,
                fillColor: "#FF5500",
                fillOpacity: 0.5,
                label : "${ad}",
                fontSize: "8px",
                fontFamily: "Courier New, monospace",
                labelXOffset: "0.5",
                labelYOffset: "0.5"
            }})
        });
        map.addLayer(vl2);

       function handler(request) {

        vl2.addFeatures(geojson_format.read(request.responseText))
        map.zoomToExtent(bounds);
        }
        var request = OpenLayers.Request.GET({
        url: "myurljson&srsName=EPSG:23036",
        params: {},
        callback: handler
    });

please show me some way


texts:
Is it possible to have data backup? If so how long?

Except Google Earth is there any other authentic system that can provide cultivations image data of a field for last few years?


texts: Aligning Polygons to Margins in Data Driven Pages

texts:
I'm creating a data driven page of polygons (forest subcompartments) on a basemap in ArcMap 10.4.1. The scale required for outputted maps is predefined (1:5,000). Output is multi page pdfs.

Users of the outputted maps/pdfs will be required to add legends to the right hand side of the maps/pdfs (pasted jpg/bmp).

In Layout view is there a way to align the polygons so the western most portion of the polygon boundary is displayed near the left margin of the layout, i.e., keep the map pages from centering on the polygon?


texts: How to insert polyline features with the list of vertices in another polyline?

texts:
I want to insert polyline features by another polyline vertices using arcpy.
I try to list X and Y coordinates of the source layer then create array from the XY and insert the array to the new feature class. the new feature class has records and length but not show any poyline.I have three segments in the source layer but the result is more than 3 records. I don't know what's the part of my code is wrong ?

# Imports
import arcpy
# Spatial reference of input feature class
SR = arcpy.SpatialReference(4326)
InFc=r"D:\_GIS\VerticesToLine.gdb\Line"
output = r"D:\_GIS\VerticesToLine.gdb\test"
for row in arcpy.da.SearchCursor(InFc, ["OID@", "SHAPE@"]):
    print("Feature {}:".format(row[0]))
    partnum = 0
    vertexArray = arcpy.Array()
    for part in row[1]:
        for pnt in part:
            if pnt:
                long = pnt.X
                lat = pnt.Y
                point = arcpy.Point(long,lat)
                print point
                vertexArray.add(point)
                with arcpy.da.InsertCursor(output, ("SHAPE@",)) as cursor:
                    polyline = arcpy.Polyline(vertexArray, SR  )
                    cursor.insertRow((polyline,))
        partnum += 1


texts: error converting shapefile to geoJSON with Ogre online tool

texts:
I am using on online tool built with ogr2ogr to convert shapefiles to geoJSON:
http://ogre.adc4gis.com/

It's worked well except for one shapefile, and I can't see why...
The shapefile has only six attribute columns: ObjectID, a shape field, a date field and three text fields..when I try to convert all that I get is:

{"errors":["Unexpected token ,"]}

I tried selecting 'skip failures' in the convertor, and I have looked for some inconsistency in the attribute names and columns but I don't see anything unusual.

What could be giving me this odd result?


texts: How to unlisten an event in OpenLayers 4?

texts:
I am trying to suppress the single click event that adds popups to the map. I have the following function:

function addFullScreenPopups(){
      map.on('singleclick', function(evt) {
        var popup = new ol.Overlay.Popup({insertFirst: false});
        map.addOverlay(popup);
        var prettyCoord = ol.coordinate.toStringHDMS(ol.proj.transform(evt.coordinate, 'EPSG:3857', 'EPSG:4326'), 2);
        console.log("Pretty Coordinates : "+prettyCoord);
        popup.show(evt.coordinate, '

Coordinates

' + prettyCoord + '

');
    });
    }

I used to use unByKey to unlisten the single click event. However, the developers removed this method in OpenLayers 4. Therefore, I am bound to use map.un method. I am trying to use this method, but it is not doing anything. I am doing  map.un("singleclick", addFullScreenPopups()); for unlistening this event, but it is not doing anything. Why?

Update: I also looked at the signature of .on method and it is given as on(type, listener, opt_this). So I changed my .on method to  map.on('singleclick', function addFullScreenPopups(evt) {//statements}); and called  map.un('singleclick',  addFullScreenPopups()); Unfortunately, this didn't help me out.


texts: Is the loop in this python code breaking because of duplicate values in the column?

texts:
This code is meant take a shapefile in the current mxd, loop through all of the values in a collumn of said shapefile. For each value it will create a definition query, zoom to features then save a copy of the mxd.  When I run the code it completes, but only generates one mxd when I'm expecting about 25.

Is the loop breaking due to duplicate values in the Field_1 collumn? If so how can I work around, if not what is breaking the loop after one single iteration?

import arcpy, os
from arcpy import env
#Overwrite features
arcpy.env.overwriteOutput = True
shape = r"C:\NewProjects\shapes\Meters.shp"
output = r"C:\NewProjects\shapes"
layer = r"C:\NewProjects\shapes\newlyr"
mxd = arcpy.mapping.MapDocument(r"C:\NewProjects\shapes\Try.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
inset = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr1 = arcpy.mapping.ListLayers(mxd, "Meters", df)[0]  # pick the layer you want to "extract" data from
#For each work order:
rows = arcpy.SearchCursor(lyr1, "", "", "FIELD_1", "FIELD_1 A")
for row in rows:
    name = str(row.getValue('FIELD_1'))
    sqlExp = """"FIELD_1" ='"""  + str(row.getValue('FIELD_1'))+"""'"""    # Create a clause to select only the current record
    lyr1.definitionQuery = sqlExp
    df.extent = lyr1.getSelectedExtent(True)
newMxd = "Figure1_" + name
mxd.saveACopy(r"C:\NewProjects\shapes\\" + newMxd + ".mxd")
print "Saving ", newMxd + ".mxd"
print "All done."
del mxd


texts: How do I teach ogr2ogr about a projection?

texts:
I have JSON generated from ArcGIS that I want to convert to GeoJSON in EPSG 4326.  This fails:

$ ogr2ogr -f GeoJSON -t_srs EPSG:4326 out.geo.json in.arcgis.json OGRGeoJSON
ERROR 6: EPSG PCS/GCS code 102100 not found in EPSG support files.  Is this a valid
EPSG coordinate system?

I found a forum post that says 102100 can be changed to 3857, which is mathematically identical.  So, this works, with the same error message:

$ ogr2ogr -f GeoJSON -s_srs EPSG:3857 -t_srs EPSG:4326 out.geo.json in.arcgis.json OGRGeoJSON
ERROR 6: EPSG PCS/GCS code 102100 not found in EPSG support files.  Is this a valid
EPSG coordinate system?

However, I have several files in the 102100 projection system.  I'd like to teach ogr2ogr about the 102100 projection, so I can script the conversion.  Is there a way to update the config so that 102100 is recognized and treated the same as 3857?

Update - This data is from the REST API of ArgGIS Server 9.3.1.  If I don't specify an output spatial reference, then I get 102100, which they define as the projected coordinate system "WGS_1984_Web_Mercator_Auxiliary_Sphere".  A much easier way to accomplish my goal is to ask for the desired projection, with the outSR=4326 query parameter.  The question is still valid if you don't have access to original data.


texts: External hyperlinks to control Openlayers map

texts:
I have a basic map defined like show below (pseudo code).

After the map.addLayers call there is a displayData(theid) which I want to call from the main body of the html document, the problem is I can't call the method.  I can achieve what I am after if I bind an onclick event to the id of the  tag but there's going to be too many links to bind to each and every  tag

Clicking on one of the hyperlinks below results in no method displayData() defined.

Thanks for reading!

View on map

View on map

View on map

View on map


texts:
I'm on my reasearch for GPS and sensors integration. I've just wondered that to integrate between GPS and camera you need time offset synchronize which called as a timestamp.

But how  does software calculate those time offset estimate automatically? (I use APM - Mission Planner for software and GPS module on arduino for hardware), I have googled it but I haven't found what I want.

Because I'm working in "accuracy and precision quality" of integrated sensors, what I want to find is like "Is there any fluence that take to the distance among the sensors as the affect of time offset?" or "let say there is, how can I calculate that in concrete way? Like if the time offset goes one second, how much the distance lack (in meters)".

Can anyone give me some explanation about this?

Or do you have any paper or article about it?

Or maybe have a right keyword for google search?


texts:
With ArcEngine 9.3.1 in java (but I suppose, it's the same in .NET), I'm developping an application which loads a mxd file. This file contains about 20 layers with shapefile, JPG's, TIFF's files. This file is located on a remote server but accessible through windows share.

After instantiating a MapBean, I try to access some information on this object but I get an Exception "Argument is not a com object :  null".

Here an extract of my code:

//some code initializing the engine

mapBean = new MapBean();
mapBean.setDocumentFilename("myMXD.mxd");
add("Center",mapBean); //adding into javax.swing.JComponent

//some Intructions
int count=mapBean.getLayerCount();

//count is 0 but after some time, it becomes the right value !!

So I tried to wait for some thread having finished their jobs in such way:

Thread thread;
Thread[] threadArray = new Thread[100];
int threadCount = Thread.enumerate(threadArray);
for (int threadIndex = 0; threadIndex < threadCount; threadIndex++)
{
    try {
        thread = threadArray[threadIndex];
        if (thread.getName().startsWith("RetrieveThread")
                || thread.getName().startsWith("Image Fetcher") )
        {
            thread.join();
            threadIndex = 0;
            threadCount=Thread.enumerate(threadArray);
       }
            } catch (Exception ex) {
                ex.printStackTrace();
            }
        }

But it doesn't really help.

Has someone an idea how to know when the MapBean is really ready with all his information ?

Thanks in advance,

Fabien


texts: Seeking st_multi equivalent in Oracle Spatial or converting polygon to multipolygon

texts:
Is there an equivalent to the PostGIS ST_MULTI available in Oracle Spatial?

or

How can I convert POLYGON to MULTIPOLYGON in Oracle spatial?


texts: Google Earth / QGIS projection (great circle line)

texts:
Here is my question:

1- Based on a long distance segment ( approx 2000 km) drawn on Google Earth and imported as shapefile in QGIS I want to obtain a terrain profile based on "Terrain Profile" plugin.

2- The problem is that my great circle segment drawn in Google Earth is not going over the same terrain in QGIS. This can be seen clearly using the Google earth layer in QGIS:

Google Earth:

QGIS:

3- The difference is approx 45 km between the 2 lines and the projection of my raster is WGS84 (EPGS:4326).
I tried to change my raster projection to 900913 (Google Mercator) but this doesn't work either.

Can you please let me know if there are any way to make my very long google earth segment coincide in QGIS ?

Many thanks in advance


texts: ogr2ogr transform failure: "illegal axis orientation combination"

texts:
I'm attempting to use ogr2ogr to reproject a MA community boundaries shapefile (available here).  Question: what part of the command causes the transformation fail? Is either my code or the source projection invalid?

Let's call the TOWNSSURVEY_POLY.shp shapefile src.shp for convenience. My command,

ogr2ogr -f "GeoJSON" -t_srs EPSG:4326 dest.json src.shp

yields the following error for which which I cannot find any documentation:

ERROR 1: illegal axis orientation combination
ERROR 1: Failed to reproject feature 0 (geometry probably out of source or destination SRS).
ERROR 1: Terminating translation prematurely after failed

src.shp has the following projection:

PROJ.4 : '+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs '

OGC WKT :
PROJCS["Lambert_Conformal_Conic",
    GEOGCS["GCS_GRS 1980(IUGG, 1980)",
        DATUM["unknown",
            SPHEROID["GRS80",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",41.71666666666667],
    PARAMETER["standard_parallel_2",42.68333333333333],
    PARAMETER["latitude_of_origin",41],
    PARAMETER["central_meridian",-71.5],
    PARAMETER["false_easting",200000],
    PARAMETER["false_northing",750000],
    UNIT["Meter",1]]

While dest.json is created, it contains no features.

Edit: yep, I meant EPSG:4326, and have changed it in the code above. The error remains despite the fix.  cat dest.json yields an empty string.  gdal-config --version yields 2.1.0 (I use the anaconda gdal dist). I will attempt to (re)compile gdal 2.1.2 from source tomorrow, as it is the most recent stable source release.
 I've edited in the exact shapefile name as well.
Final edit: completely forgot about this.  Iirc, the coordinates were flipped.


texts:
This is getting worse...
I created a project with about 6 vector layers (which I digitised) with some raster tiles.
It is stored locally on my Windows 7 PC, running QGIS 2.0

The project was fine until I came to open it yesterday, when I got the Bad Layers dialogue. Entering the correct paths does not fix it.

I've opened up various iterations of this project I had made, and my backups & they're all broken. Some I definitely experimented with Tools/Scripts (IRO my other post), but I don't think I had in the earlier ones. One of these projects' "bad layers" is a temp output from a tool.

Not sure if the time coincides, but the list of Recent Projects has reverted to that of some weeks ago.  I did a System Restore, which assured me my data was safe, but as I'm investigating, it suggests this is the culprit - if the temporary file in AppData has gone that would trip it up.

Where is the "recently used files" list - is that in AppData too? (this is less of a problem, just a possible clue to the above)

I've skimmed the .qgs file in Notepad++ to see if I could identify these layers, but can't work it out. Can this be edited like the old ArcView  .apr files?


texts:
I'd like to try to optimize a world-wide 10m contour lines (isohypses) Postgis database(from SRTM).
The complete DB takes 504GB on disk, contour lines are already in epsg:3857, and simplified with ogr2ogr at 4m. The DB is used for rendering tiles with Mapnik.
What I tried so far is:

convert elevation to smallint instead of numeric(12,3).
add a boolean column to 50m and 100m lines instead of using the % operator in mapnik queries, indexed or not.
CLUSTER using ST_GEOHASH().

The results are shown below for a 20x20km and a 50x50km box obtained on a 7200rpm HDD, and I'm a bit disappointed cause I expected huge improvements.
These test results are obtained on a small 630MB extract (~550x750km in the Alps).
I'm not surprised the clustering did not add much improvement: the contours lines have been pushed into Postgresql sequentially from 1x1° shapefiles, so the disk order is probably fine enough.

Question: Maybe I made mistakes in my method, or maybe I forgot to try something else ? Any suggestions ?
More details below.
Convert elevation from numeric() to smallint:
ALTER TABLE contours ALTER COLUMN height TYPE smallint;

Boolean column for 50m lines:
ALTER TABLE contours ADD COLUMN _50m boolean;
UPDATE contours SET _50m = NOT CAST( height % 50 AS BOOLEAN );
CREATE INDEX _50m_geom_idx ON contours (_50m);

Cluster by geohash:
DROP INDEX contours_geom_idx;
CREATE INDEX contours_geohash ON contours (ST_GeoHash(ST_Transform(wkb_geometry,4326)));
CLUSTER contours USING contours_geohash;
DROP INDEX contours_geohash;
CREATE INDEX contours_geom_idx ON contours USING GIST (wkb_geometry);

A typical query (adapted from mapnik queries):
EXPLAIN ANALYSE SELECT ST_AsBinary("wkb_geometry") AS geom
FROM
(SELECT wkb_geometry,height::integer FROM contours WHERE _50m) AS "contours-50"
WHERE "wkb_geometry" && ST_SetSRID('BOX2D(698000 5678000,748000 5728000)'::box3d, 3857);

Drop cache between each test query:

sudo sync;
sudo service postgresql stop;
sudo echo 3 | sudo tee /proc/sys/vm/drop_caches;
sudo service postgresql start


texts: GeoJSOn or WFS service to view, edit and save geometries

texts:
To view, edit and save in postgis geometries using OpenLayers + Geoserver
what is recomenadado use GeoJSON or serivços WFS?
I ask this for performance and additionally want the edition is recorded in the database realtime.


texts: QGIS is freezing on a transformation between EPSG:26912 and EPSG:3402.  It appears to be looking for a grid shift file AB_CSRS.DAC which I can get from the federal government but I don't know how or where to apply it in QGIS.  What am I missing?


texts: relating land cover raster to vector data

texts:
I have a land cover raster basemap that covers different land cover/uses in Europe to 100m resolution and a separate Excel legend file.
Separately I have a point shapefile with the details and location of archaeological sites, most of the sites fall into the landuse category of agricultural land, but I want to show this more clearly with a simplified map, and to eventually display a table showing something like
"Pastoral Land 3 sites
Arable land 10 sites" etc

I vectorised the raster at first, but am not sure if that was the right thing to do. I am playing around with habitat evaluation plugins, but what I really want to find out is a technique of evaluating the relationship between the vector point (site) and the basemap (landuse), rather then going through the rather arduous process of manually zooming in and noting the position of the point and the colour of the pixel, then going through the legend.  To make the final result clearer on a printout display Id also like to change the image resolution.

Is there a specific plugin I can use for this work? I also am trying to reclassify the raster as the detail now is more than I need.

I am using QGIS 2.8 Vienna. Below is an image of the two layers


texts: Label cutting off in data driven pages of ArcMap

texts:
I am using ArcMap 10.6.1.

I've been trying to play around with the Label Manager but can't seem to figure out a solution.  I am trying to label line features in Data Driven Pages and for those lines that are along the border of my map, the label is being shown as half in / half out.

Is there a way to make sure the labels stay completely in the map area, or even completely out?

The picture below shows what I'm looking at.


texts:
I converted a polygon in an mdx file I have to a shapefile, then converted its vertices to points, then used "add X Y feature" and choose coordinates "WSG 84", now when I open the dbf file I get coordinates like "3.17328051950, 3.07370669600" and they are really out of the country I am using (Egypt), so what is wrong here?

Points to help:
- The mxd file I am using has Egypt and some neighboring countries around it not the whole world..

I used the MID option when converting to point.


texts: Dragging one shapefile into QGIS creates three layers?

texts:
I am using QGIS 3.0.1.  Everytime I drag in a shapefile (actually, *.shp, .shx and .dbf) I get three layers, that seem to differ just in symbology color.  This happens with both shapefiles I get from outside sources and ones that I have generated in Python.  My colleagues using ArcGIS see just a single layer when they do the same operation.

This isn't a problem other than adding what looks like redundant layers - or maybe I do not understand something kind of basic here.


texts:
I would like to change the default path for saving layers ("save as" in the main layers panel).
Now the default is C:\Program Files\QGIS Lyon\bin (in which I even don't have the rights to save apparently...)

Many answers on this forum tell that it should not be like this, and that QGIS will save in the last folder used.
This is not the case in my situation.
I even saved the whole QGIS project in the by me desired folder.
Every single SHP file I want to save, is by directed to the QGIS Lyon\bin folder by default; unless I first click the "browse"-button, then it goes to the desired folder; If I don't first click the "browse"-button, it directs to the bin folder.


texts: How to get attribute table information when clicking on OpenLayers webgis?

texts: Might seems like a very trivial question to ask, but I am currently facing the problem with webGIS using JavaScript and python (with flask and psycopg2) where I load some shapefiles into pgadmin with PostGIS installed. I can see the vector layers, but I cannot think of solution how I can get information from attribute table.


texts: Putting 2 maps together under condition in Google Earth Engine?

texts:
I'm new in GEE.

I've made 2 maps: dNBR and NBRpost

and I just want to make 1 resulting risk map that put both together under the condition:

dNBR below threshould1 Not burned
NBRpost below threshould2 Moderate

NBRpost above threshould2 Extreme

I started to make filters like:

var filt1=(ee.Filter.lt(dNBR, -2));
var filt2=(ee.Filter.lt(NBRpost, -233));

Where do I go from here?


texts: Extract longitude and latitude from shapefile

texts:
I want to create a media wiki page with map option where polygon required the following format. What is name of this format and how can it possible to extract this format from shapefile?

polygons=62.103882522897855,5.09765625:58.309488840677645,5.712890625:58.95000823335702,10.8984375:61.68987220046001,12.83203125:63.35212928507874,11.865234375:64.1297836764257,13.974609375

Detailed format can be found here

http://semantic-mediawiki.org/wiki/Maps_examples/Google_Maps_shapes


texts: Microsoft SQL Server - Multi-user solution

texts:
My team and I are currently trying to apply multi-user techniques in our projects. So far we have used an ArcGIS Enterprise build for this purpose but we quickly learned there are a few limitations. (One of them being lack of topology).

As a solution to our needs, we have installed an SQL Express server with a workgroup licence. I did a test on a small project with two editors and it seemed to work fine but after that, another team complained about performance. (The server does not lack resources and it is local)

Is SQL a good practice for multi-user editing in ArcMap? If so, is it worth leaving the users to edit on the DEFAULT version rather than a replica or another version?


texts:
I have url for XYZ source tiles. But there are file names like:

/1/1-0.jpg, --> {z}/{x}-{y}.jpg

There is old code for ol2 transform url.

MyClass = OpenLayers.Class(OpenLayers.Layer.XYZ, {
getURL: function (bounds) {
    var xyz = this.getXYZ(bounds);
/... code .../
    var textY = xyz.y.toString(2);
    while (textY.length < xyz.z) textY = '0' + textY;

    for (var i = 6; i < xyz.z; i++)
            path += textY.substr(0,i-5) + '-' + textX.substr(0,i-5) + '/';

    path += textY + '-' + textX + '.jpg';

    return path;
},

CLASS_NAME: "MyClass"

It is possible to make it in ol3?


texts:
I've tried the x86 and x64 versions of OSGeo4W, but gdal2tiles is missing.

Has gdal2tiles been removed from OSGeo4W?

               -={ OSGeo4W Shell Commands }=-

   bmp2tiff            msencrypt           scalebar
   cs2cs               nad2bin             shp2img
   curl                nccopy              shptree
   dllupdate           ncdump              shptreetst
   fax2ps              ncgen               shptreevis
   fax2tiff            ncgen3              sortshp
   gdaladdo            nearblack           testepsg
   gdalbuildvrt        nircmd              textreplace
   gdaldem             nircmdc             thumbnail
   gdalenhance         ogdi_import         tiff2bw
   gdalinfo            ogdi_info           tiff2pdf
   gdallocationinfo    ogr2ogr             tiff2ps
   gdalmanage          ogrinfo             tiff2rgba
   gdalserver          ogrlineref          tiffcmp
   gdalsrsinfo         ogrtindex           tiffcp
   gdaltindex          opj_compress        tiffcrop
   gdaltransform       opj_decompress      tiffdither
   gdalwarp            opj_dump            tiffdump
   gdal_contour        osgeo4w-setup       tiffinfo
   gdal_grid           pal2rgb             tiffmedian
   gdal_rasterize      portmap             tiffset
   gdal_translate      ppm2tiff            tiffsplit
   geod                proj                tile4ms
   geotifcp            psql                vcredist_2013_x86
   gif2tiff            ras2tiff            xmlcatalog
   gltpd               raw2tiff            xmllint
   legend              rgb2ycbcr           xxmklink
   listgeo             rpcgen
   mapserv             rpcinfo

   apache-install     o-help             setup
   apache-restart     o4w_env
   apache-uninstall   setup-test

GDAL 2.1.0, released 2016/04/25


texts: Tables not showing within schemas following Table retrieval in QGIS (Oracle Spatial)

texts:
When connecting to an Oracle Spatial database, QGIS is not consistently populating the list of tables under each schema. When I try to expand a schema (as I have tried to do with PARKING in the picture below), nothing shows.

If I can expand a schema whilst the table retrieval is occurring then sometimes a list of tables will populate and I can access their contents but tables are missing from the list.


texts: How to call 'shapely' .coords on a fiona filtered list?

texts:
I have to filter a shapefile by attribute then compare it to another shapefile using .coords...however I am a bit perplexed how to do this as my filter is returning a list type...(and .coords requires a shapely Geometric object).

Is there another way to filter my shapefile and then get the coordinates of the filtered items?

import fiona
from shapely import Point
# open a shapefile of points
with fiona.open("file.shp") as src:
   # returns a list - w all attributes and point geoms
   filter_list = filter(lambda f: f['properties']['foo']=='bar', src)
   # need to run .coords on filter_list


texts:
I'm looking for a way to customize the layout of popup.
What is asket to me is to change the defined labels of navigation inserting a custom string and move and change the icons to navigate

the first change is something like this:
change the header title
"(1 of 2])"
to
"element 1 of 2"

the second change is to move the navigation tools from the header to the footer and to use labels instead of images (something like previous/next)

Is that possible or do i need to implement a own popuptemplate?
thanks


texts: Copying Grass data from one location to another in Grass GIS

texts:
I have created several Grass locations with single mapsets while processing a series of satellite imageries at different stages of my analysis.

Now I need to copy some preprocessed rasters and vectors from the older locations to my current location and mapset. I am aware of using 'g.copy' which can be used in case of a single grass location and different mapsets. I know I can also export the rasters and vectors and import into my current location.

But found it wasting time as those data are already imported in Grass and all these Grass Locations are in same 'region' settings and with same projections.

Any hack for adding those rasters and vectors in my current Location and mapset?

USing Grass GIS 7.4 on Windows 7.


texts: Plot animation between source/destination IP addresses

texts: Can anyone share code with me that will plot and animate "attack lines or points" between a source IP address and destination IP address, or vice versa.  I'm trying to recreate the Norse dashboard here: http://www.norse-corp.com/technology_overview.html


texts: Get indices from pixels with a certain value in PyQGIS

texts: I'm working in the Python-Console of QGIS and I have an object form the class QgsRasterLayer. I want to get the indices of those pixels that have certain pixel value (e.g. I need the indices of all the pixels with a pixel value = 50), but I don't know how.


texts: Inferring lat/lon based on known location coordinates

texts:
I have a point dataset with a set of coordinates that is in an unknown projection. For some of these points I know the unprojected latitude and longitude. Based on these known data points, would it be possible to infer what the latitude and longitudes are for the unknown points? Is there a formula that can do this?

Here is a sample of the data - I would like to get the Lat/Lon of the first five based on what I know for the last 5 data points:

PROJECTED_X PROJECTED_Y Latitude        Longitude
2635536     1326596
2253740     2804378
2561993     1430129
2919984     1356273
2919931     1355933
2658206     1371867     33.8105204      -117.8349355
2363687     1478454     34.11562481     -118.8053659
2101778     2290419     36.35547224     -119.6548331
2095540     2274891     36.31321472     -119.6764159
2699433     1470349     34.07804599     -117.691484


texts: Create paper topography anaglyph map that can be visualized in 3D using red-green glasses

texts: I've recently seen paper maps showing elevation displayed using a blue-red gradient that could be visualized in 3D with cheap paper red-green glasses. The effect is quite neat and I'd like to learn how to add this same effect to my maps using QGIS. Is it possible?


texts: GDAL set negative values to nodata or NaN

texts:
How can I set all negative values in a raster to NaN or no data?

I've made a virtual raster to link data and coordinates from a HDF5 then using something like this:

gdalwarp -geoloc -tr 3000 3000 -srcnodata -9999 -dstnodata -9999 -r average -s_srs EPSG:4326 -t_srs '+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m' input1.vrt output1.tif


texts:
I am using QGIS.

I have one quite large map (A0 or A1) and I want to print in like atlas (just split it to print on A4 or A3). If understand correctly, I have to use Atlas functionality in Composer.  I have done all steps from official manual (inc. creation of coverage grid), but by not very clear reason I am getting 81 almost equal pages. On each page region is just slowly moving up. What am I doing wrong?

Ok, I want to clarify what I want to do. I have one huge map and I want to split it into submaps to print on small papersheets (and I want to specify its size by hand). Each submap should have its own frame with grid values, legend and so on. Or suppose you have detailed map of Earth and want print it like ordinary paper atlas.  I think this can be done without working on each separate page. But I can't figure out how to do this in QGIS. But am not an expert:), I am just trying to create my first map for backpack trip.
If this feature is not implemented, but nice to have, I can try to write a plugin if I have time.

I have never used ArcGIS, but I found this http://help.arcgis.com/EN/ARCGISDESKTOP/10.0/HELP/index.html#//00s90000002s000000.htm . Looks quite close to my idea.


texts: Aggregating polygons with Basic level licence of ArcGIS Desktop?

texts:
I need to aggregate near-by polygons, but it has to work with a Basic ArcGIS Desktop Licence. I can't use the Aggregate tool (only Advanced) and I can't use the buffer tool, because the rounded edges create a lot of unnecessary points (flat edges are not available with a Basic License).

The accuracy is not very important (the less points created the better), it should look something like this:


texts: Pythonaddins Message Box encoding problem

texts:
I've been trying to create a MessageBox with accented characters.  To do so, I used the decode() and encode() functions. The code is part of a Python Toolbox.

pythonaddins.MessageBox("ééàà".decode('utf-8').encode('cp1252'), "My title", 0)

And that's the result I get:

I've also set the default encoding to utf-8 at the beginning of my script:

`# -- coding: utf-8 --

Strangely, when I try it directly in the ArcGIS Python console, it gives the right result:

It seems like running the code from the Python console in ArcGIS doesn't give the same result as running it from a Python tool ...

I've tried lots of things to fix this like:

Setting the text string in Unicode

pythonaddins.MessageBox(u'ééàà'.encode('cp1252'), "My title", 0)

Using Windows popups

ctypes.windll.user32.MessageBoxA(0, "ééàà", "My title", 0)

... and it still behaves the same way (works in ArcGIS Python console but not from the toolbox).

Ideas why this is happening?

EDIT:
This is an example of tool code that doesn't give the right result:

# -*- coding: utf-8 -*-
import arcpy
import pythonaddins

class MsgBoxTest(object):

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

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

class MsgBoxTest(object):

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

def getParameterInfo(self):
    """Define parameter definitions"""
    params = None
    return params

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

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

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

def execute(self, parameters, messages):
    """The source code of the tool."""
    message = "ààéé"
    pythonaddins.MessageBox(message.decode('utf-8').encode('cp1252'),"Attention",0)
    return


texts:
I installed the Customize ToolBars plugin and created a few custom toolbars. However, I saw they only get a label name instead of the actual toolbar name. I checked customization and also printed all toolbars via Python:

from PyQt4.QtGui import QToolBar

for x in iface.mainWindow().findChildren(QToolBar):
    if x:
        print x.objectName()

This results in the following outputs. The plugins I want to find names for are showing up blank or empty.

Is there a way to just add some name to these unnamed toolbars?


texts:
Here is what I am doing:

Create a layer by uploading a shapefile into QGIS
Select attributes of that file
Export selected attributes with the "save only selected features" check box when I save the file

I need to export names corresponding to each polygon. Unfortunately when I save as KML file and open in Google My Maps, I can see only polygon boundaries, therefore I need to name them. The info is in a column in the attribute table.


texts:
As total newbie of QGIS, I am slowly progressing into using properly QGIS but I lack of knowledge of the software.
I am missing the knowledge to edit .shp files in the goal of adding or changing data.
Context: I have a shape layer already being done and I would like to add data in each specific areas so I can fetch later the information of the area with the 'Identity features' tool.
Finally, last question is there an efficient way to drop the data from an Excel file (.xls) to the QGIS shape layer(s) ? The context would be that I have some updated data and after making a new Excel with all of the updated data for each parcel I would drop them into the QGIS to update the map so I can simply click on the desired zone with the 'Identity features' tool to have the needed data.
This is an example of a selected area of the shape layer

And this is the tab I would like to be able to edit

I added some pictures to make sure my rambling is more clear because I was using QGIS in French before so some terms might be different (or poorly translated by myself).
I've been also trying to follow guides on the matter but it's a bit complicated with the mix of outdated ones and the up to date.


texts:
I'm trying to join a field named "Area" belonging to the shapefile 'Ancient Woodland' with the shapefile 'Net Development Area' (see screenshots below). The results only return 'NULL' which I don't understand because the field type of "Area" is a number. Even if I join the entire table to Development area all the entries return 'NULL'

Both shapefiles have a common field which I'm using to perform the join ("Substations"). Surely I don't need to create a .csvt to do this?


texts:
I have an app that uses slippy maps.  The tiles I'm currently using contain a lot of information that is superfluous to my users.

I've been investigating generating replacement map tiles that contain only very basic features such as country/province boundaries, rivers, lakes, EEZ lines (the application focuses more on water than land) and maybe railways and major motorways when at high zoom levels.

First I tried loading planet.osm.bz2 (a 21GB file) into Maperitive.  This appeared to be a bad idea as the process memory usage rose quickly beyond that which my machine could tolerate.  I've tried again on smaller files (such as Antarctica alone) and have more success, but I really need to generate tiles for the entire Earth.

Most of the features in these large OSM data files are not of interest to me.

Can someone please recommend an approach that should run on my x64 Win7 PC with 4GB RAM?


texts: Filter by day and reducing region by polygon using image collection in GEE

texts:
My overall program is here. I have an aerosol index that I already reduced for daily means:

var countries = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017");
var Iraq = countries.filter(ee.Filter.eq('country_na', 'Iraq'));

var airq_19 = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_AER_AI')
  .select('absorbing_aerosol_index').
  map(function(image){return image.clip(Iraq)})
  .filterDate('2019-01-01', '2019-01-31');

var daysOfMonth = ee.List.sequence(1, 31); //Restricted to first 31 for trial run

var composites = ee.ImageCollection.fromImages(daysOfMonth.map(function(d) {
  var filtered = airq_19.filter(ee.Filter.calendarRange({
    start: d,
    field: 'day_of_year'
  }));
  return filtered.mean().set('day', d);
}));

print("Composite aerosol", composites);

I also have a grid over Iraq:

//Bounding box
var lat_start = 29.0990251735;
var lat_end = 37.3852635768;
var lon_start = 38.7923405291;
var lon_end = 48.5679712258;

//Set number of grid cells
var num_cells = 10;
var lon_edge = (lon_end-lon_start)/Math.sqrt(num_cells);
var lat_edge = (lat_end-lat_start)/Math.sqrt(num_cells);
// 3) Create the grid
var polys = [];
var polys_line = [];
var cell_id = 0;
for (var lon = lon_start; lon < lon_end; lon += lon_edge) {
  var x1 = lon;
  var x2 = lon + lon_edge;
  for (var lat = lat_start; lat < lat_end; lat += lat_edge) {
    cell_id = cell_id + 1;
    var y1 = lat;
    var y2 = lat + lat_edge;

    polys.push(ee.Feature(ee.Geometry.Rectangle(x1, y1, x2, y2), {label: cell_id}));
  }
}

//Final grid
var grid = ee.FeatureCollection(polys).filterBounds(Iraq);
print('Final grid', grid);
Map.addLayer(grid);

I am trying to build a function that identifies mean aerosol by day within each grid. I have developed a function that is currently not running. It looks like this:

var daysOfYear = ee.List.sequence(1, 365);
var counts = grid.map(function(feature) {
  var counts = daysOfYear.map(function(day) {
    var dayEvents = composites
      .filter(ee.Filter.eq("day", day))
      .filterBounds(feature.geometry())
      .reduceRegions ({
    collection: grid,
    reducer: ee.Reducer.mean(),
    scale: 100
      });
    var key = ee.String('day').cat(ee.String(ee.Number(day).int()))
    return ee.List([key, dayEvents.size()]);
  }).flatten()
  return feature.set(ee.Dictionary(counts))
});

print("aerosol by grid and day", counts.limit(10));

The function above results in the error:

composites.filter(...).filterBounds(...).reduceRegions is not a function

Is the issue that the aerosol composites are stored as an image collection? If no, then how can I modify the function above such that the function will produce the intended feature collection?


texts: Creating multiple polygon features based on the number of point features contained

texts:
Software: ArcGIS Desktop 10.x or ArcGIS Pro

I have a polygon dataset, and a point dataset. Any number of points (or zero points) may intersect with the polygon. I need to generate a new polygon dataset which contains polygon with the attributes of points which intersect them. Ordinarily this would just be a spatial join, however I need a separate polygon feature for each intersecting point.

I need this to be achieved via a new dataset, rather than relating the point table to the polygon table.

I feel like there is a built-in way to do this in ArcMap, but can't place it.


texts:
I am using a QgsScaleComboBox and it's behavior is really strange. I also have 2 buttons in my gui for adding/deleting items. Those are the methods I use for adding/deleting:
Adding:

    def addScale(self):
        newScale = (str(1. / self.dialogui.scaleBox.scale()))
        if self.dialogui.scaleBox.findText("1:" + newScale) == -1:
            self.dialogui.scaleBox.addItem("1:" + newScale)

Deleting:

    def deleteScale(self):
        delitem = self.dialogui.scaleBox.currentIndex()
        self.dialogui.scaleBox.removeItem(delitem)

I also store all items in the QgsScaleComboBox after the plugin closes. I store them like this:

    list = []
    for i in range(self.dialogui.scaleBox.count()):
        list.append(self.dialogui.scaleBox.itemText(i))
    QSettings().setValue("scales", ";".join(list))

And when the plugin starts I get the values back like this:

        if not QSettings.value("scales") is None:
            for scale in QSettings.value("scales").split(";"):
                self.retrieveScale(scale)

And with this method I add them back to the QgsScaleComboBox:

    def retrieveScale(self, scale):
        if self.dialogui.scaleBox.findText(scale) == -1:
           self.dialogui.scaleBox.addItem(scale)

I also have some predefined scales that I add with the updateScales method, it looks like this:

self.dialogui.scaleBox.updateScales(predefinedScales)

Now if I add a new item it adds the item to the QgsScaleComboBox but I can't select it. If I print the value of the current selected item in QgsScaleComboBox with the currentText() method it prints the right item. What am I missing?


texts:
I had developed a plugin that counts the number of features for a layer in QGIS 2.6. I am using Windows 7 os. It works fine when i loaded a layer. But when I remove that added layer and then if I click the plugin it gets crashed. The snapshot is attached.

The part of the dump file error is also attached:

.

Here is function which does the main work of the plugin.

def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()
    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.featurecombo.addItem( layer.name(), layer )
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result  == 1:
        index = self.dlg.featurecombo.currentIndex()
        layer = self.dlg.featurecombo.itemData(index)
        if layer > 0:
            QMessageBox.information(self.iface.mainWindow(),"feature count","%s has %d features." %(layer.name(),layer.featureCount()))
        else:
            QMessageBox.information(self.iface.mainWindow(),"Warning!",'Touch me after loading a layer!')
            # Do something useful here - delete the line containing pass and
            # substitute with your code.
            pass


texts: Adding field in attribute table on several layers using ArcGIS ModelBuilder?

texts:
I have about twenty separate polygon-layers, all consisting of almost - but not entirely, the same information in their current attribute tabels. However, now I want to add a field in the attribute table all of the layers.

I don't want to merge the layers and I'm thinking that there must be a smoother way to do it beside adding one field at a time in each layer.

Can I perhaps use the ModelBuilder?

Update:

Right now, I'm just testing with one layer, thinking I can use the same model later. But for some reason this doesn't turn out quite right.

The first steps works fine, adding fields to the original layer. After that I'm intersecting the layer with three other layers to be able to calculate the fields, but this is were something goes wrong and it's like the field calculator calculates the wrong fields for some reason.


texts: Missing Scottish Administrative Data in OS Opendata sets

texts:
I'm not sure if this forum is really the correct place to ask, but impertinence would lead me to believe that there may be some one on here that may know the answer.
I'm using the Ordnance Survey Open Data bdline_gb shape files to try to pull together some geographic coverage data for my organisation that covers England Scotland and Wales.
I'm working mainly with the Parish data set; because that's how ridiculously specific my organisation can be (I am determinedly prising them away from part parishes!). And fortunately England and Wales administrative Boundaries "build-up" from parishes (or communities in Wales).
In the bdline_gb data set there is also very nicely the District, Borough, County and Unitary Authority shape files as well; so I can overlay to see where I have coterminous boundaries that no longer need to be down to Parish level.

Unfortunately Scotland doesn't [still] have Parishes; but additionally the data sets don't seem to include any ADMINISTRATIVE boundaries for Scotland, at least not in it's proper coded sense; just the filler polygons without the appropriate codes.
I can get electoral data till it's coming out of every orifice just not Administrative-
does anyone know where I can source this?