texts: Unifying color scale of different tiles in Google Earth Engine?

texts:
I found some other posts asked the similar questions but in different tools like ArcGIS, GDAL, QGIS etc. What I want to ask is how to do this in GEE.

Now I have removed the bad observation pixels (like clouds and scan-line corrector-off gaps in southern Florida. But I found every tile has its own unique color scale. Particularly for the tile Path:16 Row:43 (the one locates at the southwestern corner which covers Key West) has a significant difference with other tiles, which results in a big error when I calculating and filtering the indices.

Here I attached the code and screenshot for your reference.

//Choose country using GEE Feature Collection

var region = ee.Geometry.Polygon({
  coords: [[[-82.48, 26.93], [-79.91, 26.93], [-79.97, 25.68],
            [-80.78, 24.75], [-81.21, 25.07]]],
  geodesic: false
});

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;

var non_mosaic = landsat8 .filterBounds(region) .filterDate('2017-01-01', '2017-12-31');

var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// A function to mask out cloudy pixels.
var cloud_shadows = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
  // Return an image masking out cloudy areas.
};

// A function to mask out cloudy pixels.
var clouds = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 5,5, 'Cloud').eq(0);
  // Return an image masking out cloudy areas.
};

var maskClouds = function(image) {
  var cs = cloud_shadows(image);
  var c = clouds(image);
  image = image.updateMask(cs);
  return image.updateMask(c);
};

var mosaic_free = non_mosaic.map(maskClouds).median();

var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};

Map.addLayer(mosaic_free, visParams, 'Cloud free');


texts:
I am using the following code and wanted to know how I can change it from click event to hover or mousemove event. I know there is a eventListner but I am unable to locate the click event where I must change it to mousemove event.

info = new OpenLayers.Control.WMSGetFeatureInfo({
            url: 'http://localhost:8080/geoserver/iPlant/wms',
            title: 'Identify features by clicking',
            queryVisible: true,
            eventListeners: {
                getfeatureinfo: function(event) {
                    map.addPopup(new OpenLayers.Popup.FramedCloud(
                        "chicken",
                        map.getLonLatFromPixel(event.xy),
                        null,
                        event.text,
                        null,
                        true
                    ));
                }
            }
        });
        map.addControl(info);
        info.activate();

Update 2.0

I made the changes and also added the PopUp, but only the click function is being recognized.

info = {
            click: new OpenLayers.Control.WMSGetFeatureInfo({
            url : 'http://niles.iplantcollaborative.org:8080/geoserver/wms',
            title : 'Identify features by clicking',
            queryVisible : true,
            eventListeners : {
            getfeatureinfo : function(event) {
            deletePopUp();
            popup = new OpenLayers.Popup.FramedCloud("popInfo", map.getLonLatFromPixel(event.xy), null, event.text, null, true);
            map.addPopup(popup);}}
            }),
            hover: new OpenLayers.Control.WMSGetFeatureInfo({
            url : 'http://niles.iplantcollaborative.org:8080/geoserver/wms',
            title : 'Identify features by clicking',
            queryVisible : true,
                    hover : true,
            eventListeners : {
            getfeatureinfo : function(event) {
            deletePopUp();
            popup = new OpenLayers.Popup.FramedCloud("popInfo", map.getLonLatFromPixel(event.xy), null, event.text, null, true);
            map.addPopup(popup);
            }}
            })
            };
            map.addControl(info.hover);
            info.hover.activate();


texts: Arcpy Spatial Join Fails: Field is Not Nullable

texts:
We are working to perform a spatial join of a very large point .shp of voters (target_features) with another of incidents (join_features).

arcpy.SpatialJoin_analysis(target_features, join_features, out_features, "JOIN_ONE_TO_MANY")

is working on a small subset of the voters, but when run on the entire set, the error is:

Field is not nullable

for each of the fields in the incidents .shp.  Because this is Arc 10.2, we cannot AlterField_management or change the Nullability in Properties as seems possible in 10.3.  Instead, we copied the .shp features to a Default.gbd Feature Class (which now claims to have nullable fields in Properties) and re-ran it, but received the same error.

Is there another way to change or ignore the nullability of fields?


texts: Shapely LineString and Polygon intersect?

texts:
here is my python code, it does not intersect and when I plot this in google map , I see the intersection. shapely.geometry is in Cartesian plane, how you are solving this?

l = "40.65:-114.07 40.21:-112.96 40:-112.66 39.4:-112.7 39.39:-113.28 39.68:-113.94 40.42:-114.26 40.6:-114.22"
    points = [[float(x) for x in c.split(":")] for c in l.split(" ")]
    print(points)
    polygon = Polygon(points)
    path = LineString([(42.049999,-96.25), (38.0880979565,-119.609216457)])
    if path.intersects(polygon):
        print("intersect")

here is my code for google map:

https://gist.github.com/anonymous/8393151


texts: Updating row in list of feature classes using dictionary in Python?

texts:
I have a list of feature classes (fcList) that I would like to do the following.
 - add a field (fire_dis)
 - update row from a single value from a dictionary(fcnaDict)
 - more specifically I am trying to run this through a loop so that the first item in the fcList will get the first dictionary value "Livingston Fire" and so on until every fc has the reciprocal dictionary value.

Here is my code.  FYI I am quite new at python.

a = arcpy.MakeFeatureLayer_management(outFeatureClass, "LivingstonFire", CLF)
b = arcpy.MakeFeatureLayer_management(outFeatureClass, "GardinerFD_GatewayHoseCo", GFDGHC)
c = arcpy.MakeFeatureLayer_management(outFeatureClass, "CookeCitySilverGateVFD", CCSVFD)
d = arcpy.MakeFeatureLayer_management(outFeatureClass, "ClydeParkCity", CPC)
e = arcpy.MakeFeatureLayer_management(outFeatureClass, "ClydeParkRFD", CPRFD)
f = arcpy.MakeFeatureLayer_management(outFeatureClass, "WilsallRFD", WRFD)
g = arcpy.MakeFeatureLayer_management(outFeatureClass, "ParkCountyRFD1", PCRFD1)
h = arcpy.MakeFeatureLayer_management(outFeatureClass, "ParadiseValleyFSA", PVFSA)
i = arcpy.MakeFeatureLayer_management(outFeatureClass, "CountyProtection", CoPro)
j = arcpy.MakeFeatureLayer_management(outFeatureClass, "NullParcels", Null)

a2 = arcpy.Dissolve_management(a, str(a), "", "", "MULTI_PART","")
b2 = arcpy.Dissolve_management(b, str(b), "", "", "MULTI_PART","")
c2 = arcpy.Dissolve_management(c, str(c), "", "", "MULTI_PART","")
d2 = arcpy.Dissolve_management(d, str(d), "", "", "MULTI_PART","")
e2 = arcpy.Dissolve_management(e, str(e), "", "", "MULTI_PART","")
f2 = arcpy.Dissolve_management(f, str(f), "", "", "MULTI_PART","")
g2 = arcpy.Dissolve_management(g, str(g), "", "", "MULTI_PART","")
h2 = arcpy.Dissolve_management(h, str(h), "", "", "MULTI_PART","")
i2 = arcpy.Dissolve_management(i, str(i), "", "", "MULTI_PART","")
j2 = arcpy.Dissolve_management(j, str(j), "", "", "MULTI_PART","")
#Dissolve the exempt layers from above
k2 = arcpy.Dissolve_management(OtherExlayer, "OtherExparcels","", "", "MULTI_PART","")
l2 = arcpy.Dissolve_management(FSlayer, "FSparcels","", "", "MULTI_PART","")
m2 = arcpy.Dissolve_management(DNRClayer, "DNRCparcels","", "", "MULTI_PART","")

## This is my dictionary
fcnaDict = { "a2":"Livingston Fire", "b2":"Gardiner Fire Department & Gateway Hose Co",
"c2":"Cooke City Silver Gate VFD", "d2":"Clyde Park City", "e2":"Clyde Park RFD",
"f2":"Wilsall RFD", "g2":"Park County RFD1", "h2":"Paradise Valley FSA",
"i2":"County Protection", "j2":"Null Parcels", "k2":"Other Exparcels",
"l2":"FS parcels", "m2":"DNRC parcels"}

## This is my list of objects for feature classes
fcList = [ a2, b2, c2, d2, e2, f2, g2, h2, i2, j2, k2, l2, m2]

## This is my loop with the above mentioned items
for fc in fcList:
    arcpy.AddField_management( fc,"fire_dis", "TEXT", 100)
    cursor = arcpy.da.UpdateCursor(fc,"fire_dis",)
    for row in cursor:
        row[0] = fcnaDict.values()
        cursor.updateRow(row)

The error I am getting now is that updateRow does not support lists.


texts:
Here is the code block (the layers are stored in Oracle database):

var saveStrategy = new OpenLayers.Strategy.Save();

var adresspt_wfs = new OpenLayers.Layer.Vector("Address points WFS", {
     strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
        protocol: new OpenLayers.Protocol.WFS({
        url: "http://localhost:8080/geoserver/wfs",
        featureType: "RAILWAY",
        featureNS: "http://localhost:8080/MIA",
        srsName: "EPSG:900913",
        geometryName: "SDO_GEOMETRY"
      })
    });

map.addLayer(adresspt_wfs);

Please correct me if I'm wrong: It adds new vector layer "Address points WFS", on which we can edit/delete/add some vector objects. After save all that will be saved RAILWAY table as a geometry objects in SDO_GEOMETRY column using WFS service, that is part of GeoServer (in my case I use GeoServer)... Is that right?

Also, I don't quite understand what is "featureNS" and where to find the proper one for my layer.

If I open: http://localhost:8080/geoserver/wfs/DescribeFeatureType?version=1.1.0&typename=MIA:RAILWAY  in address bar, the browser hangs on opening that schema...

MIA is my default workspace, its URI is http://localhost:8080/MIA

I've tried many configurations for initializing the WFS protocol, but can't save drawn geometry into database... Please show me what am I doing wrong during this initialization.

Please consider that I'm completely new to "geo" programming.


texts: How to fix the reprojection from EASE-2 grid product SMAP to geographic coordinates?

texts:
I have been working with SMAP satellite data, especially for moisture and soil properties.

I follow the idea of use GDAL solve everything, and make something similar to this published in Link to first approach to download smap data

Modifing the code and testing:

import os
import h5py
import numpy as np
from osgeo import gdal, gdal_array, osr

# the file to download from   server

https://n5eil01u.ecs.nsidc.org/SMAP/SPL4SMAU.003/2017.08.01/SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5

path = "/path/to/data"
h5File = h5py.File(path + "SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5", 'r')
data = h5File.get('Analysis_Data/sm_rootzone_analysis')
lat = h5File.get("cell_lat")
lon = h5File.get("cell_lon")

np_data = np.array(data)
np_lat = np.array(lat)
np_lon = np.array(lon)

num_cols = float(np_data.shape[1])
num_rows = float(np_data.shape[0])

xmin = np_lon.min()
xmax = np_lon.max()
ymin = np_lat.min()
ymax = np_lat.max()
xres = (xmax - xmin) / num_cols
yres = (ymax - ymin) / num_rows

nrows, ncols = np_data.shape
xres = (xmax - xmin) / float(ncols)
yres = (ymax - ymin) / float(nrows)
geotransform = (xmin, xres, 0, ymax, 0, -yres)

dataFileOutput = path + "sm_rootzone_analysis.tif"
output_raster = gdal.GetDriverByName('GTiff').Create(dataFileOutput, ncols, nrows, 1, gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

output_raster.SetProjection(srs.ExportToWkt())
output_raster.GetRasterBand(1).WriteArray(np_data)  # Writes my array to the raster

del output_raster

So, using this approach, the result is a global map with many problems of projections, as for example the image below, produced by the python code above.

To compare with a correct data, the same image was extract from h5, using HEG nasa software.

first suggestion

Following the modification suggested by @Ali:

The solution

Proposed by @AndreJ in Solution that solve the problem


texts: Union large number of rasters in PostGIS

texts:
I'm having some difficulty finding the sum of a large number of rasters in PostGIS. Basically I have nearly 1 million polygons which I converted to rasters in a PostGIS table. I now want to overlay all of them and create a single output raster which is the sum of all their pixel weights.

I have tried using something like

SELECT ST_Union(rast, 'SUM') FROM public.my_raster_table

This seemed to work to work well on a small data set, but the pixel values in the output were not what I expected. Many hundreds of these rasters overlap, and so I expected the range of pixels to be around 0-500. Instead they only ranged from 1-6. In other words, if seems that if you try to sum more than a handful of rasters with this method, data is lost.

When I tried the method on the full set of data, it ran extremely fast (<2 seconds), which seems impossible for adding nearly 1 million rasters. The output appeared to be completely flat; i.e., pixel values ranged from NaN to 1. What is going on?

Update: Apparently some of the rasters I was using were unexpectedly huge (e.g., entire western half of United states at a pixel size of ~1/3 mi. It is clear that these are what was causing the server disconnect, though it's still unclear why. I'm also not sure what to do about it, except to filter out the really gigantic ones.


texts: QGIS 3 map viewing area not filling up window area

texts: I just installed QGIS 3.2.1 from QGIS 2.18 on MacOS 10.3.6. When I open the new install, the map does not fill up the area in the app window (below).


texts: Moving row down or up manually in attribute table of ArcMap?

texts:
I have digitized several lines in ArcMap and I did not pay any attention to ordering it initially only to realize later that I need to other them somehow not according to any of the columns.

Is there a way to manually re-order the rows?

I have tried doing it in editor mode and outside it but did not work.


texts:
There are many ways to preview a map as it is seen by people with color vision impairment. ColorOracle, preview mode in QGIS, etc.
But none of these options allow exporting the modified image, only printscreen and save.
How can I transform map into color vision impairment sight?


texts: How to filter date values using OGR 1.11.0?

texts:
With OGR 1.10.1 this works:

ogrinfo test2.shp -sql 'SELECT COUNT(*) FROM test2 WHERE DATE < "1970-01-01"'

Note that DATE is a column of type date (not string).

(Here's test data if you want to try: http://ge.tt/6DBclJh/v/0)

The same query in 1.11.0 causes

OGR[3] error 1: Type mismatch or improper type of arguments to < operator.

Has the syntax changed or is this really a bug?

I filed a ticket, so in case this is a bug, I'm all set: http://trac.osgeo.org/gdal/ticket/5507


texts: Create S-57 chart on ARCGIS containing only shapefiles

texts:
I have a requirement where a shapefile/files needs to be converted to S-57 ENC chart for use on ECDIS system. I am not working in GIS domain so dont necessary understand all technical bits and bobs involved.

Till now I have done following steps.
I have downloaded ARCGIS + Maritime Charting and Production mapping (Desktop 10.4.1) on trial basis. Please note that I dont have any ENC data to work on.

I have created a new file Geodatabase in Acrmap, Imported XML workspace data/schema from file "Nautical_chart_template_GX.xml".

I have carried out  steps to load the shapefiles under "UserdefinedfeaturesA", create AOI (using Feature builder) under "PL_ENC_AOI" setup scale to 10000000.

Start "New Nautical Product" wizard, which works without any errors.  But when creating ENC "Export Nautical product". The I am getting "Error - Nothing found 10000000" during Export process. Export log does not give much details.

I dont have any data inside other layers. Just PL_ENC_AOI to define AOI and "UserdefinedfeaturesA" with 3 sublayers contains data.

Can anyone guide me the right way to approach this issue.


texts: Reload attribute table does not work for attribute changed by dataProvider

texts:
In QGIS2.18, I'm changing feature attributes using the following two methods:

dataProvider:

attrs = { 0 : "hello", 1 : 123 }
layer.dataProvider().changeAttributeValues({ fid : attrs })

layer object:

layer.startEditing()
layer.changeAttributeValue(fid, fieldIndex, value)
layer.commitChanges()

I keep the attribute table open, and monitor the changes. The first method dataProvider is not able to update the attribute table on the fly, even the "reload table" button on the top menu of the table does not work. I have to re-open the table to see the changes. The second method is working, the attribute table gets instantly updated.

I would like to stay with the dataProvider method but not sure how to get instantly updated attribute table?


texts: How can VIIRS-Data be obtained / processed

texts:
I am looking for a way to obtain data about avg light-pollution for specific places.
It seemed NASAs https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A1/ might be a good source.
However I am totally lost on their websites. Does anybody know any starting-points / documentation on:

How to download the data in any way.
How to extract the images for regions from that data. Is that what https://blackmarble.gsfc.nasa.gov/Tools.html does?


texts: How to quickly add a DWG with multiple layers?

texts:
In QGIS 2.14 "Essen", I need to add a DWG with multiple layers.

So I choose "add vector layer" feature in QGIS 2.0 to import the .dwg/.dxf file.

However, QGIS adds all layers from the DWG file into one layer in QGIS.
How to import the DWG so each layer will be in a separate QGIS-layer;
That is, each layer in QGIS should be named the same as layer is called in the DWG file.


texts: Creating empty list and adding items to this newly created list in Earth Engine Code Editor?

texts:
I have tried something like this:

var empty_list = ee.List([])

However, when adding items to this list, it returns an empty list regardless.

Here is a code snippet. Note that "col_list" is a variable representing an image collection.

var list = ee.List([]);

for(var i = 0; i < 5; i++){

  var img = ee.Image(col_list.get(i));

  var dictionary = ee.Dictionary({
    coords: img.geometry().getInfo()["coordinates"],
    landsat_id: img.get('LANDSAT_ID'),
    satellite: img.get('SATELLITE'),
    sensing_time: img.get('SENSING_TIME'),
    system_index: img.get('system:index')
  });

  var feature = ee.Feature(null, dictionary);
  var list = list.add(feature);
}

This definitely looks wrong but I am not sure how else to do so.


texts: I have sucessfully installed AniMove by activating experimental plugins (yay! thanks gis.stackexchange). But I can't find it in Qgis! Yes I have activated it. Where is it? Is this a common bug? I searched for it the toolbox, I looked up the names of the several single tools based on this youtube video and only found "Kernel Density", which is filed under SAGA. But I don't have SAGA installed (and don't know if I need to?)


texts: Define custom CRS in WKT from point and angle

texts:
I'm trying to define a custom CRS using the WKT Syntax. However when I do the projection I'm off by about 2km.

Here is my rotation point.

Local X and Y:
X: 4635.396 Y: 2397.085

MGA94 Zone50:
x: 560255.527
y: 7427753.462

Control Points:

Mine X | Mine Y| MGA 94(50) X| MGA 94(50) Y
2453.122|3210.002|563053.406|7431461.771
-1735.225|-853.24|557798.872|7428929.256
5663.648|7386.58|567416.171|7434410.368
12607.859|-1438.839|571218.306|7423848.605
8502.84|2620.24|568605.287|7428993.832
-2500.032|3457.767|558433.331|7433259.449

These are the steps I'm following base on WKT for local mine grid:

Convert the MGA94 Zone50 (EPSG:28350) x and y to longitude and Latitude ("EPSG:4326"). I've used the python package pyproj

from pyproj import Transformer
transformer = Transformer.from_crs("EPSG:28350", "EPSG:4326", always_xy=True)
print(transformer.transform(564420.896, 7430150.547))

This gives the points
(117.62970383981178, -23.236582623614485)

Base on the link above I've put used the WKT synta

PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",
GEOGCS["GCS_GRS 1980(IUGG, 1980)",
DATUM["D_unknown",
SPHEROID["GRS80",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
PARAMETER["latitude_of_center",-23.25839260829391]
PARAMETER["longitude_of_center", 117.58908484003899],
PARAMETER["azimuth",-18.39841101],
PARAMETER["scale_factor",0.999585495],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1]]

I then paste the code into QGIS custom CRS.

When I apply this custom CRS to a polygon layer in QGIS the polygon appears about 2km away from the actual location.

Can anyone offer any advice on how to achieve more accuracy?


texts: Dividing an irregular polygon into sub-regions

texts:
I am trying to divide an irregular polygon into equal widths and am having trouble understanding the proper tool to use. I am in possession of ArcGIS and am aware that the Parcel Editor is the closest function to accomplish this task, but I am not actually able to edit said polygon - any ideas? I am also aware that I could divide the regions manually with the cut polygons tool, but I was looking for something more automatic.

Below is my goal, mainly to divide the regions into approximately an equal width of w.


texts: Getting maximum values of attribute of multiple Shapefiles merged into new shapefile in QGIS?

texts:
I have thousands of points in multiple shapefiles. I want to get the maximum value of an attribute of each shapefile. Then I want to merge all the maximum values to a new layer.

Is there a way to extract max values as a batch process from multiple layers?

EDIT:

This is the sample attribute table. The max values I need are from "DTM50austri" from each Shapefile.

This is the query I've tried. "Zusammengeführt" is the merged shapefile.


texts: Find if a table is mappable/contains an Obj column

texts:
How do I check whether a table is mappable?

I am trying to select all records from a table that contain objects, like this:

Select * From TableName Where Obj Into TempTable

If a table doesn't have an Obj column then this will cause an error. How can I check for the existance of the Obj column?


texts: Optimising performance on local QGIS + PostGIS installation?

texts:
I'm learning how to manage a local PostGIS (ver 2.1) database, with spatial indexing enabled on the shapefiles that I've imported (using the PostGIS Shapefile Import/Export Manager). I have imported a large shapefile with 4.7 million polygons and a 2.1GB attribute table into a schema.

However, using QGIS 2.8.2 to load, render, and open the attribute table of this table has been rather disappointing. Adding the dataset to the QGIS canvas freezes the application for a short period. While the polygons are quickly rendered, opening the attribute table produces a small dialogue window saying '____ features loaded' until all 4.7 million features are loaded. While QGIS takes a long time to open the attribute table, it also appears to consume a great deal of memory (~1.1GB) as well. I've avoided making queries to subset rows of the attribute table due to poor performance as well.

I've heard of other QGIS+PostGIS users reporting excellent performance with tables containing >500,000 or even millions features, but maybe my expectations are unreasonable. So while I'd like to enjoy the performance benefits of PostGIS, I'm not sure what the problem is. Here's what I've tried:

I've created a spatial index on the table using the using the PostGIS Shapefile Import/Export Manager and the following SQL code:

CREATE INDEX parcels_geom_idx
   ON parcels
   USING gist
   (geom);
In QGIS under the Data Sources tab, I have the 'Attribute table row cache' set to 10,000 (the default)
The database is on the local host machine which runs on a Core 2 Duo, 8GB RAM, and a 240GB SSD

Eventually I'd like to store up to 100GB of spatial data from a collection of geodatabases in a PostGIS database, but with performance problems on just a few data sets I'm not confident in moving away from ESRI's platform towards an open source ecosystem of GIS software.


texts: Matching paths based on both GPS pathways and the respective time recorded

texts:
I have two database tables containing individual lat-long records of vessels moving around in waters for a given time period:

Database A contains a set of records of a particular type of vessel that engages in a range of activities
Database B contains a set of records of all types of vessels moving around in a particular area of interest

The objective of this analysis is to find the distribution of the time taken in undertaking a certain activity for the vessels in database A.

I have to make a few operational assumptions that take place in reality for this particular activity, but in essence, the matching that I would like to execute between these two datasets would be based on the two key criteria between two vessels in database A and B:

Coinciding GPS paths
Same time frame

I have attempted to look through at a few other map matching questions and answers but the proposed solutions appear to revolve around methods that do not take into consideration the factor of time,  (i.e. I require the paths of two vessels to coincide at the same time) and that the GPS paths needed to be matched solely on coinciding paths (with some degree of inaccuracy to be allowed in the GPS coordinates).

The ideal output of this analysis would be a list of vessels from database A that has "matched" with a vessel from database B for a particular time frame (starting time and ending time as two columns), the starting and ending lat-long locations (the coordinates where their paths begin coinciding, and the coordinates their paths diverge).

The recording frequency of each vessel's records depends on the speed the vessel is traveling in the waters.  Higher travel speeds would result in a higher frequency of records.  This means that the recording does not take place in regular intervals, and the time would need to be matched within specific time frames.

As I am really new to using GIS in solving the above issues, I might not be aware of all the correct terminology I should be using in referencing my data or any specific methodology in problem solving in this field.  Please excuse me, and let me know if any details have been missed out, or if anything is not clear enough.

I am platform agnostic in finding a solution to these, so please do point me in the right direction should there be something already available in QGIS/Python etc.


texts:
I'd like to turn Google Earth Engine geometries to feature collection, however receive error message like:
Collection: Geometries cannot have their properties modified or be placed into collections.

How should I modify the code?
My code is following:
var glc = ee.FeatureCollection('users/wangq/test_slope')
var geometries = glc.geometry().geometries()
print(geometries)

var shape = ee.FeatureCollection(geometries)
print(shape)


texts: Zoom out grayed out in ArcGIS API for JavaScript?

texts:
I'm using ArcGIS API for JavaScript 3.13

The zoom out slide button is grayed out, however it works when clicked. Does anybody know how to make it not grayed out?

My map initialization code looks like this:

esri.config.defaults.map.zoomDuration = 0;
esri.config.defaults.map.zoomRate = 0;
_map = new Map("myMap", {
     center: [-100, 48],
     zoom: 3,
     basemap: "streets",
     minZoom: 3,
     navigationMode: "classic",
     logo: false
});

jsfiddle here:
https://jsfiddle.net/xnqt6z37/7/

I figured out that it is the   "minZoom" attribute that's making the zoomOut button grayed out. minZoom is set to be the same as initial zoom, so initially map cannot be zoomed out further.

Also,

esri.config.defaults.map.zoomDuration = 0;
esri.config.defaults.map.zoomRate = 0;

is causing the problem, I added these to remove zoom animation due to a bug in chrome, and zooming with  ESRI arcgis:
Zoom hangs ArcGIS API for JavaScript application?

I cannot find a setMinZoom() method in map, even though there is a getMinZoom().

I have this as a hacky workaround:
https://jsfiddle.net/xnqt6z37/13/
it sets the initial zoom to be smaller than minZoom, then set zoom back to the desired initial zoom, which is the same as minZoom in this case.

Any better ways?


texts: Are US Census tracts GEOID the same across different years?

texts:
I want to know if it is appropriate to join tables from different years to the same census  tracts TIGER file.

If all census tracts' IDs do not stay the same, then how would they change when adding new tracts, modification, merging, etc?


texts: Generate GeoJSON tile layers from Shapefile

texts:
I have some vector layers in a shapefile that I would like to display in the Google Maps Javascript API. I recently learned that you can generate GeoJSON from shapefiles in the correct CRS using QGIS: Convert Shapefiles to GeoJSON? that worked well for me, but the problem is that the file is large - the resulting JSON is 250MB in my case. So I'd like to have it as a collection of tiles instead. I discovered that the Google Maps API can load tiles from GeoJSON as described here
https://stackoverflow.com/questions/27284251/vector-tiles-with-google-maps-v3
But I can't find a way to actually generate tiled GeoJSON like that. I noticed that both QGIS and MapTiler have an option to write vector tiles in the XYZ structure, but the individual files come out as .pbf, which Google Maps can't handle. So I am looking to create them as GeoJSON. Is there a way to do that?
Edit: In searching around, I also discovered tippecanoe and was able to run it, but it looks like it also generates pbf files. I want the individual tiles in geoJSON.
Edit 2: I then discovered that I could run tippecanoe-decode on the very same pbf files that tippecanoe itself had created. While a little weird, it seemed to work and generated JSON that looked like below, which is throwing errors in Google, I'm assuming because of the nested FeatureCollection...
    {
  "type": "FeatureCollection",
  "properties": {
    "zoom": 18,
    "x": 42016,
    "y": 91728
  },
  "features": [
    {
      "type": "FeatureCollection",
      "properties": {
        "layer": "sea",
        "version": 2,
        "extent": 4096
      },
      "features": [
        {
          "type": "Feature",
          "properties": {
            "autocad_el": 0,
            "autocad_th": 0,
            "autocad_wi": 0,
            "autocad_co": 106,
            "autocad_la": "C_1CON_ROAD_CURB",
            "autocad_li": "ByLayer"
          },
          "geometry": {
            "type": "LineString",
            "coordinates": [
              [
                -122.299832,
                47.442148
              ],
              [
                -122.299805,
                47.44213
              ],
              [
                -122.299647,
                47.442021
              ],
              [
                -122.299621,
                47.442003
              ]
            ]
          }
        },
        {
          "type": "Feature",
          "properties": {
            "autocad_el": 0,
            "autocad_th": 0,
            "autocad_wi": 0,
            "autocad_co": 106,
            "autocad_la": "C_1CON_ROAD_CURB",
            "autocad_li": "ByLayer"
          },
          "geometry": {
            "type": "LineString",
            "coordinates": [
              [
                -122.299832,
                47.442153
              ],
              [
                -122.299805,
                47.442135
              ],
              [
                -122.29964,
                47.442021
              ],
              [
                -122.299613,
                47.442003
              ]
            ]
          }
        }
      ]
    }
  ]
}


texts: Managing Layer's Visibility in DynamicMapServiceLayer iOS

texts:
I'm working on a iOS mobile application using Arcgis SDK.
There are number of layers(and sublayers) in my Dynamic map service.I want to turn ON/OFF visibility of selected layers(as well as their sublayers) depending on user's choice.
I tried doing so by setting 'visibleLayers' property of my AGSDynamicMapServiceLayer.But it works only once i.e. when I turn ON any layer's visibility, it gets visible (along with some other layers/sub layers as well). After that when I try to turn OFF this layer's visibility,It doesn't do so.


texts: How can I align the pixels with using GDAL and Python?

texts:
I am trying to use Python/GDAL to align different GeoTIFFs that I warped (gdal.Warp) to the same SRS. My problem is one of the images is rotated, the pixels do not overlap perfectly, which is something I need.
How can I do this? I am fine with deleting the outermost pixels for resampling, I just do not know how to do it. Here is an example, as you can see the orange one is slightly off (~1.5m in x and y at 10m resolution)

Here is my current code:
    Da = gdal.Open(i)
    Name = i[7:]
    options = gdal.WarpOptions(srcNodata=-1000,dstNodata=-1000,dstSRS="EPSG:32632")
    DB = gdal.Warp(destNameOrDestDS = "Band02/rep_cut/"+Name[:-4] + "_32632.tif",srcDSOrSrcDSTab= Da, options = options)

If add targetAlignedPixels = True to the options, there is no output at all.


texts: Setting up Python and QGIS for standalone applications

texts:
I am learning QGIS for an upcoming internship. I have been trying to configure Python 3.7.5 with QGIS 3.8 (a Coruna) for two weeks now. I have overcome quite some issues by reading a lot of posts, but now I seem to be stuck.

I am trying to run a simple testScript.py from a windows 10 terminal and got many many errors about modules not being found etc. Solved those issues by reading a lot of posts here and elsewhere, but now I finally get access to my script's imports and also managed to get a Python GUI to open up, the cmd complains about "ERROR 4: Unable to open EPSG support file gcs.csv"

Allthough,  as you can see, I did follow the hint and I do point to the folder where this file can be found:

Microsoft Windows [Version 10.0.17134.1069]
(c) 2018 Microsoft Corporation. All rights reserved.

C:\Users\Kunga Led>SET OSGEO4W_ROOT=C:\OSGeo4W64

C:\Users\Kunga Led>SET QGIS_PREFIX=%OSGEO4W_ROOT%\apps\qgis

C:\Users\Kunga Led>SET PATH=%PATH%,%QGIS_PREFIX%\bin

C:\Users\Kunga Led>SET PYTHONHOME=C:\OSGeo4W64\apps\Python37

C:\Users\Kunga Led>SET PYTHONPATH=%QGIS_PREFIX%\python;%PYTHONPATH%

C:\Users\Kunga Led>SET QT_QPA_PLATFORM_PLUGIN_PATH=C:\OSGeo4W64\apps\Qt5\plugins

C:\Users\Kunga Led>SET GDAL_DATA=C:\OSGeo4W64\share\epsg_csv

C:\Users\Kunga Led>cd Desktop

C:\Users\Kunga Led\Desktop>python testScript.py
ERROR 4: Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

How do I solve this issue?

I have read that the OSGEO4W64 installer actually takes care of these issues by itself, but since I installed it, I haven't had anything but trouble getting stuff to run. (Of course I am totally new to this so there will be a lot I need to learn about, I just don't know exactly where to look) I would like to add that running scripts from the built-in QGIS console or editor is never a problem. The problems occur when I try to run a .py file from the commandline.


texts: How can I use osm2po to calculate distance matrix in Km? There is an example that calculates distance matrix in time but I cannot make it turn into kilometers. I dont understand how should I use the int[] path parameter into the calcPathLength function.


texts: Is there a On Map Move event in OpenLayers?

texts: In my OpenLayers web app I have registered an onmousemove event which performs certain functions.  What I am really after though is a way to register and handle the user panning or zooming the map.


texts: Processing time/speed for Mosaic to New Raster

texts:
I have a question on speeding up “Mosaic to New Raster” or finding an alternative.  After completing a suitability analysis we are trying to define movement corridors.  We used a modified CircuitScape that left 4 folders, each containing 1080 raster tiles across our study area.  We now need combine all 1080 rasters in each folder then combine these 4 combined rasters into a final raster which will be analyzed for paths.

I tried a python script that just wouldn’t work so I went on to Model Builder.  It iterates through the grids in the folder, collects the values (raster names), then mosaics all to a new raster.  It has been running for 5 days now and the results window indicates it is working.  I’m hoping it will conclude in the next day or 2.  But, we still have 3 more folders, each with 1080 rasters, that need to be mosaicked.  So….is there a way to speed up Mosaic to New Raster or is there an alternative I could use to produce the final grid?

The cell size is 270 meters, 540 x 540 columns and rows, 1 band, 32 bit floating.


texts: Finding tie points between a georeferenced and a non geo referenced image

texts: Is there an open source program for finding tie points between a non-georeferenced image and a georeferenced image? I am looking for something I  can run at the command line in linux to incorporate into a processing chain, that outputs the tie points in image coordinates in the un-referenced image,  as well as in geographic coordinates.


texts: ogrmerge.py - How to running from Python

texts:
I need merge list of shapefiles.
It is possible using ogrmerge.py
I use it and run using subprocess:
import pathlib
import subprocess

class MergeShpFiles:
    def __init__(self, output_file: pathlib.Path):
        self.output_file = output_file

    def merge_shp_files(self, file_list: list):
        command = "ogrmerge.py -single -o "
        command += str(self.output_file) + " "
        for file in file_list:
            command += str(file) + " "
        "-src_layer_field_name country"
        subprocess.call(command, shell=True)

It is possible to run function in Python without subprocess?
It is problem when I switch between Windows and Linux


texts: Cannot create layer because "unable to open database file" in QGIS 3.2.0?

texts:
I'm brand new to QGIS but 10+ years ago I was adept with ArcGIS. I installed 3.2.0 just a few days ago and all went well at first. I located/downloaded/imported the data I want and added a field to the table in each layer. I'm using Windows 10 Home v.1803, 64-bit. Now I want to create new layers using that new field. Screeching halt.

I can't create a new layer. I can't clip the layer. I can't filter the layer. I can't create anything else, apparently, because I keep getting some version of this error:

  "Export to vector file failed. Error: Creation of data source failed
  (OGR error: sqlite3_open(NC GeoRegions.gpkg) failed: unable to open
  database file)"

I'm not sure what database it's trying to open, or I would check out its read-write permission. And anyway I can't believe the message is accurate. If QGIS cannot "open" a database, how in the world can it read that database to populate attribute tables? It makes no sense on its face.

What am I doing wrong?


texts: Finding GRASS module after upgrading QGIS (Mac)?

texts:
I upgraded from QGIS 2.10.1 PISA to QGIS 2.18 last week (via KyngChaos), and since then have not been able to run GRASS within QGIS - it appears not to be able to find the module.

I had a look at this related question, and I have located the processing options (screengrab below) and tried to find the correct grass folder paths for my system, however, I can just see the QGIS application in my applications folder - no further folder structure to locate the correct file?

Amending the current file paths to the suggested paths I can find (/Applications/QGIS.app/Contents/MacOS/grass or /Applications/GRASS/GRASS-7.x.app/Contents/MacOS) results in a 'Wrong value for parameter "GRASS folder" message.

Can anyone offer any advice? I am a QGIS newbie running on Mac OS Sierra.


texts: Adding elevation to shapefile in ArcGIS Desktop

texts:
I have a shapefile that's actually a polyline (but it's 200 miles!)
on the other hand I have a SRTM layer with elevation data.

How should I add the elevation data from SRTM layer to the shapefile?

I cannot use 3D "interpolate line" and draw a line because it's 200 miles!

I am using ArcGIS Desktop 10.7.


texts: Can anybody give me some articles or show me some link where GIS is used to analyse and manage micro-environment, such as a college campus or a small locality. I need to do a project on micro-environment analysis using GIS, so i need some references from where i can get some ideas. Please help.


texts:
I would like to add in a layer's field a default value. I mean, every time a new feature is created, that field would be automatically filled with the default value.

In my case the default value would be the @project_filename variable (project's file name).

I can't find this functionality anywhere.

Regards,


texts: Proper way to rbind SpatialPolygonsDataFrames with identical polygon IDs?

texts:
What's the proper R idiom for rbinding SPDFs together when the IDs overlap?  Note that here (as is often the case) the IDs are basically meaningless so it's pretty annoying that I can't just make rbind ignore them....

library(sp)
library(UScensus2000)
library(UScensus2000tract)

data(state) # for state names
states <- gsub( " ", "_", tolower(state.name) )
datanames <- paste(states,"tract", sep=".")
data( list=datanames )
lst <- lapply(datanames,get)

nation <- do.call( rbind, lst )
Error in validObject(res) :
  invalid class “SpatialPolygons” object: non-unique Polygons ID slot values

# This non-exported function designed to solve this doesn't seem to work any more.
d <- sp:::makeUniqueIDs( list(arizona.tract,delaware.tract) )
Error in slot(i, "ID") :
  no slot of name "ID" for this object of class "SpatialPolygonsDataFrame"


texts:
I want to plot some geometrical shapes but I dont know what values to put in a vector so they can be plotted. How do you know what numbers should be inserted in a vector so they can create certain shapes?

Can you provide a few shapes such circle, triangle, etc?

More importantly is to explain how to be able to know what numbers to use.


texts: Load style with 'importNamedStyle' in QGIS3

texts:
I am saving my layer-styles in the same spatialite-database, where the layer-data is saved and I am able to load the layer (with the default style) to the canvas (using PyQGIS in a plugin).

But how do I load an alternative style for this layer from my spatialite-database (from the table 'layer_styles')?

I have found the method loadNamedStyleFromDatabase(self, db: str, uri: str), but I have troubles to set the parameters.
I think db is like G:/path/to/database.db. But how to set the parameter uri correct to load an alternative layer-style.

Or is there an other way to do this. E.g. I can get the qml- or sld-definition from the table 'layer_styles' for the wanted layer. But how to load the layer with this definitions?

EDIT on 16.05.2019:

I have now tried to load the style with importNamedStyle like this:

layer = QgsVectorLayer(uri, base_name, provider)

style_qml = layer.getStyleFromDatabase(str(style_id))
style_doc = QDomDocument("qgis")
style_doc.setContent(str(style_qml))
layer.importNamedStyle(style_doc)

QgsProject.instance().addMapLayer(layer)
self.canvas.refresh()

I get the values for uri, base_name, provider and style_id from the database and I have controlled them, they are right (incl. style_qml).

I have e.g. two styles for one layer stored in the database. On running the code above, then only one of them is loading. No one of the two is set as the default style.

How can I load the style correct?


texts: Using CRS as parameter in QGIS Graphical Modeler

texts:
The question is loosely related to Using Calculator in QGIS Graphical Modeler since it is the same model I'm talking about. I have an algorithm ("Create Grid", in that case, but it might be any other) that needs a CRS as an input parameter, which I want to be set by the use of the model.

.

But from the list of available model inputs I suspect that this is currently not supported.

I tried to use "Vector Grid", which does not need the CRS, and assign the projection later, but I run into the same problem. Also tried to extract the CRS from another layer, but I just can't connect the extracted CRS with the algorithm where it is needed.

To be a bit more specific: I want the CRS as a model parameter rather than setting it within the model (which would then mean I have to edit the model every time I need another CRS)


texts: QGIS Cloud display in browser requests API Key

texts:
I have successfully created a QGIS project, added a OpenCyclemap layer and uploaded the data to QGIS Cloud, and published the map. When viewing the map in a browser, the OpenCyclemap background has diagonal text across it saying "API Key required"
See http://qgiscloud.com/smoky/Derby_tube_map_cloud for example.

This is the same situation on more than one map. I have another map using OpenStreetmap (created a lot earlier) which doesn't have this "API Key Required" text. This might be due to the different background or the age of the map.

Previously, I don't think the OpenCyclemap based maps had the text (although can't absolutely check).

What API Key is needed?

Where can I obtain the key?

How do I "install" the key into QGIS or QGIS Cloud (or elsewhere)?


texts: .Tif to .GeoTiff Scaling issue using GDAL

texts:
I was trying to georeference a .tif file generated from a matplotlib output using GDAL's python bindings and it seems to be off-scale. My objective was to create a geotiff file from the output of a contourf function. Code samples are pasted below.

fig = plt.figure()
plt.axis('off')
plt.contour(land_reference,levels=[-9999],origin='upper', cmap='gray',extent=(77.53,81.05,28.65,31.49))
plt.contourf(Z, levels=levels,origin='upper', cmap='RdYlGn',extent=(77.53,81.05,28.65,31.49))

fig.canvas.draw()

plt.savefig(r'C:\Users\MRG16-10417\Desktop\New\testspecies.tif',transparent='true',bbox_inches='tight')

Then I use GDAL as follows,

from osgeo import gdal, osr
import numpy as np

src_filename ='C:/Users/MRG16-10417/Desktop/New/testspecies.tif'
dst_filename = 'C:/Users/MRG16-10417/Desktop/New/newspecies.tif'

#The row and column values of the tif file was used here
nrows=252
ncols=386

src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

xmin = 77.5000000002060005
xmax = 81.0833333335421997
ymin = 28.6666666666176049
ymax = 31.4999999999532037

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)

# Set location
dst_ds.SetGeoTransform(geotransform)

# Get raster projection
epsg = 4326
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()

# Set projection
dst_ds.SetProjection(dest_wkt)

# Close files
dst_ds = None
src_ds = None

The resultant image on QGIS,
As you can see the image is slightly off scale and isn't exactly overlapping on the reference .shp file of the same geographic extents and projection.

What could be the reason for this?

Is it because the source image that is incorrect, as matplotlib could be scaling down the image?


texts:
Is there any api,pulgin for generating SLD?

Note:
I want to generate SLD programatically. I have property stored in database from which i would like to generate sld. I also want to generate various different tag according to geometry type like polygon, line or point.
Currently I am going xml way.

Edit:

I am looking for all the ways that I could generate SLD. So I have tagged the question with java and openlayers. It may be client side by openlayers or server side by java. I don't want to concatenate string or genrate xml. Is there any other way then this?

I have layers with multiple columns or I have multiple tables with multiple columns. Now I want to set predefined parameters like point size point color or polygon size and color in database and create SLD.

I want to use SLD with openlayers. Generating xml is fine but I wanted to know if there are other good ways.

var sld = '<StyledLayerDescriptor version ="1.0.0" ...> more xml here </StyledLayerDescriptor>';

wms = new OpenLayers.Layer.WMS.Post(
   "name",
   "http://localhost:8080/geoserver/wms",
   {
      'layers': 'myNs:layername',
      format: 'image/jpeg',
      sld_body: sld
   },
   {
      unsupportedBrowsers: []
   }
);


texts: Using projected coordinate system and transformation gives ERROR 000151?

texts:
I'm a beginner with GIS and I'm currently working on ArcGIS 10.3.
I have two distinct shapefiles: one within the projected coordinate system "NTF_Paris_Lambert_Sud_France" and the second one within "RGF_1993_Lambert_93".

When I opened the second shapefile I did the transformation "RGF_1993_To_NTF_NTv2 + NTF_Paris_To_NTF" which allowed me to correctly superimpose both of my shapefiles. Now I want to use the "Project" tool in order to change the projected coordinate system of the second shapefile on a permanent basis (in order to proceed to analysis with tools...). That's where I have a problem: the output shapefile is displaced by around 300 meters up to the North and a little bit to the West too.

Could someone help me with the method?

This is a screenshot of the "project" tool parameters. I did not filled the "geographic transformation" line because none of the proposed transformations is working and I always receive the message with Error 000151 (saying that the geotransformations are not valid).


texts: QGIS 2.14.3: I'm having problems when exporting my map to jpg

texts:
I am using QGIS 2.14.3 and I'm having problems when exporting my project to jpg (also when trying to export to pdf...). I have several shapefiles opened, they are basically polygons on a shapefile with the boundaries of South America and a hill shade (see pic1). I can see all the layers perfectly in my Print composer window, and I have no problems exporting to image, but the final jpg. (or pdf) file, just shows a half of the polygons layers (pic2). I'm using 50% transparency for the polygons.
Can anyone help me with this?

ps: apparently this happens when exporting with a resolution of 300 dpi -or more- but no for lower resolutions! is that somehow possible?!


texts: Forest change tutorial results don't make sense compared to area of the country

texts:
I've been following this Earth Engine tutorial on calculating forest loss by country, and running that same code for multiple other countries. What I had matched the expected output that the tutorial references (ie. 348,036,295 square meters of forest lost in Congo in 2012).

https://developers.google.com/earth-engine/tutorial_forest_03

I was curious what tree cover lost each year was as a percentage of total land area in each country, but those results don't make any sense. For Congo, there are years when 20-50% of the total country area show forest loss. When I tried this with other countries, the results are worse--with forest loss sometimes over 100% of the country's land area.

I'm getting the country's land area from CIA World Factbook: https://www.cia.gov/library/publications/the-world-factbook/geos/cg.html

Obviously I'm overlooking something very basic?

My code:

var country = ee.FeatureCollection("USDOS/LSIB/2013")
        .filterMetadata('cc', 'equals', 'CF');

// Get the forest loss image.
var hansen = ee.Image('UMD/hansen/global_forest_change_2018_v1_6')

// separating out the bands
var treecover2000 = hansen.select('treecover2000');
var loss = clipped.select('loss');
var lossAreaImage = loss.multiply(ee.Image.pixelArea());

// PRINT YEARLY LOSS BY COUNTRY (SMALL)
var lossyear = hansen.select(['lossyear'])

var lossByYear = lossAreaImage.addBands(lossyear).reduceRegion({
  reducer: ee.Reducer.sum().group({
    groupField: 1
    }),
  geometry: country,
  scale: 30,
  maxPixels: 1e13
});
print(lossByYear);

var statsFormatted = ee.List(lossByYear.get('groups'))
  .map(function(el) {
    var d = ee.Dictionary(el);
    return [ee.Number(d.get('group')).format("20%02d"), d.get('sum')];
  });
var statsDictionary = ee.Dictionary(statsFormatted.flatten());

print(statsDictionary);


texts: Join feature classes from two ArcSDE geodatabases with same schema?

texts:
I am working with two ArcSDE geodatabases -- one for test, one for production. They both have the same schema. I recalculated an attribute in the test environment and wanted to push that to production. My instinct was to join the two feature classes and use the Field Calculator, but there was no way to distinguish which column belonged to which database. So I appended _TEST to the aliases of the test class. When I ran the join again, it automatically appended _TEST to the production aliases as well. [grunt] Somehow I figured out which was which and was ready to use Field Calculator but... the option was grayed out in the right-click menu for the column. All columns, in fact. So my main question is:

How is it possible to join feature classes from two ArcSDE geodatabases and copy attribute values from one to the other?

Secondly, why on earth were my production field aliases being suffixed with _TEST? It seems like ArcMap has no way of managing a join between feature classes with the same field names. Is this true?


texts: Case Insensitive PropertyName Evaluation in SLD

texts:
I have a table layer_styles created by QGIS and, inside this table, the column stylesld.

The problem is that the PropertyNames (column names) stored inside this column and the generated SLD that geoserver explorer creates (and sends to geoserver) are UPPERCASE too, but in postgis, the column names of tables are LOWERCASE.

Is there a way to transform the LOWERCASE/UPPERCASE in column name(PropertyName) inside SLD or use a case insensitive evaluation on column name (PropertyName) ?

For example:

<ogc:Filter xmlns:ogc=""http://www.opengis.net/ogc"">
  <ogc:PropertyIsEqualTo>
    <ogc:PropertyName>
      <CaseInsensitive>TENSAO_OPE</CaseInsensitive>
    </ogc:PropertyName>
    <ogc:Literal>138 KV</ogc:Literal>
  </ogc:PropertyIsEqualTo>
</ogc:Filter>

Note that CaseInsensitive tag is only an example to explain what I need.


texts: Should GDAL be set to produce GeoTIFF files with compression? Which algorithm should be used?

texts:
I have a folder of GIS data that consists mainly of GeoTIFF files.  The whole set weighs in at about 1.2 GB. I noticed that if I pack the contents into a tarball, it smashes down to about 82 MB.  I would like to check the set into a revision control system sot it can be worked on by other people and it looks like there is some space that can be squeezed out.

The GDAL GeoTIFF driver page lists plenty of options that may be used to create compressed GeoTIFF files.  There are also plenty of options that affect the way each algorithm works.

The help page does a good job at describing the options but doesn't elaborate on how to select an algorithm or the tradeoffs that are associated with the varying level of compression.  This leads to the following questions:

The pros of using compression are a dramatic savings in space.  What are the cons?  Is information lost when the image is compressed?
How should one go about choosing an algorithm and compression level.  Do some types of images lend themselves to a certain algorithm?


texts: GRASS gis: treating certain cells as starting points in mapcalc

texts:
I am attempting to calculate the vertical height from horizontal of every point on a slope. I want to use stream beds (the lowest points in my DEM) as my starting spot, with the horizontal plane originating from there. From there, figuring out the vertical height of different points upslope should be as simple as the expression:

tan(slope.map)*horizontal distance

My horizontal distance is simply the resolution of my DEM and the slope map that was derived from it (10m).

However, I need to make sure that r.mapcalc knows that I want to use the stream beds as the origin points from which height is to be measured in that expression. I have already created a stream map using r.streams, and was now hoping someone could help me figure out how to tell r.mapcalc to use that map as the origin in figuring out vertical heights.


texts: NullReference exception when casting a layer into IGxFile Interface

texts:
I'm doing a simple test of a tool that uses a dockable window in an ArcGIS Desktop Add-In. I have a button that opens the window which will be the form that interacts with my other classes. I'm trying to access the maps and layers of an mxd but when debugging I get a NullReferenceException saying my layer object is not set to an instance of an object. Anyone know why this is happening?

namespace DockableExportTopoTool
{
    public class LaunchButton : ESRI.ArcGIS.Desktop.AddIns.Button
    {
        protected override void OnClick()
        {
            UID dockableWinUID = new UID();
            dockableWinUID.Value = ThisAddIn.IDs.DockableWindow1;
            IDockableWindow dockableWindow1 = ArcMap.DockableWindowManager.GetDockableWindow(dockableWinUID);
            dockableWindow1.Show(true);
        }
    }
    public partial class DockableWindow1 : UserControl
    {
        Site site;
        public DockableWindow1(object hook)
        {
            InitializeComponent();
            this.Hook = hook;
        }
        private void btnExport_Click(object sender, EventArgs e)
        {
            site = new Site();
        }
    }
    class Site
    {
        private IMxDocument mxdoc;
        private const string quadsLayerFilePath = @"G:\USGS 24k Map Sheets.lyr";
        public Site()
        {
            mxdoc = ArcMap.Application.Document as IMxDocument;
            map = mxdoc.FocusMap;
            SetLayerReferences();
        }
        private void SetLayerReferences()
        {
            for (int i = 0; i < map.LayerCount; i++)
            {
                IGxFile gxFile = map.Layer[i] as IGxFile;
                if (gxFile.Path == quadsLayerFilePath)
                {
                    MessageBox.Show(gxFile.Path);
                    break;
                }
            }
        }
    }
}

The gxFile object is raising the exception, right now I'm just displaying a MessageBox to show whether I've accessed the layer I need.

ArcGIS 10.1 SP1 C#.NET 4.0


texts:
result_polys = sol.vector.polygon.georegister_px_df(mask2poly,
                                               affine_obj=tile_tfm,
                                               crs= 'epsg:4326')

print(result_polys.crs)

epsg:EPSG:4326

result_polys.to_file('result_polys.geojson', driver='GeoJSON')

  CRSError: Invalid input to create CRS: epsg:EPSG:4326


texts: Error ORA-00920: invalid relational operator

texts:
I am running this SQL in Oracle Spatial in postgis worked.
In the Oracle error

[Err] ORA-00920: invalid relational operator

I need to select all geometries that touch on my streets layer, ie the corners.

SELECT DISTINCT
    (b.cdidelog) AS Cruzamento
FROM
    "GIVVI_EIXO_LOGR" A,
    (
        SELECT
            *
        FROM
            "GIVVI_EIXO_LOGR"
    ) b
WHERE
    ST_Touches (A.geom, b.geom)
AND A.cdlog = 7775026
AND A.nmidelog <> b.nmidelog


texts: Identifying subgraphs within planet_osm_line table?

texts:
I want to identify subgraphs within the planet_osm_line table of the osm2pgsql-schema. those subgraphs are classified by certain tags i want to check. My first attempt was to do this recursively with a sql-function like this:

CREATE FUNCTION private.find_subgraphs() RETURNS void AS $$
DECLARE
    good_way RECORD;
    id_rec RECORD;
    next_subgraph_rec RECORD;
    id BIGINT;
    next_subgraph_id BIGINT;
BEGIN
    SELECT osm_id FROM private.good_ways WHERE subgraph_id IS NULL LIMIT 1 INTO id_rec;
    id := id_rec.osm_id;
    WHILE (id IS NOT NULL) LOOP
        SELECT nextval('private.subgraph_seq') AS subgraph_id INTO next_subgraph_rec;
        next_subgraph_id := next_subgraph_rec.subgraph_id;

        WITH RECURSIVE waysForSubgraph(osm_id, previous_way, visited_ways, done) AS (
            SELECT osm_id, public.planet_osm_line.way, ARRAY[ROW(osm_id, public.planet_osm_line.way)], false
                FROM public.planet_osm_line
                WHERE osm_id = id
            UNION ALL
            SELECT
                    public.planet_osm_line.osm_id,
                    public.planet_osm_line.way,
                    visited_ways || ROW(public.planet_osm_line.osm_id, way),
                    ROW(public.planet_osm_line.osm_id, public.planet_osm_line.way) = ANY (visited_ways)
            FROM waysForSubgraph, public.planet_osm_line JOIN private.good_ways
                ON public.planet_osm_line.osm_id = private.good_ways.osm_id
            WHERE ST_Intersects(previous_way, public.planet_osm_line.way) AND NOT done
        )
        UPDATE private.good_ways
            SET subgraph_id = next_subgraph_id
            WHERE osm_id IN (SELECT osm_id FROM waysForSubgraph);

        SELECT osm_id FROM private.good_ways WHERE subgraph_id IS NULL LIMIT 1 INTO id_rec;
        id := id_rec.osm_id;
    END LOOP;
    RAISE NOTICE 'subgraph building done';
END;
$$ LANGUAGE plpgsql;

I know, that the assignment of id and next_subgraph_id is not perfect, but works.

Now to the problem:
For a small number of ways (like 10-20) this seems to work just fine. But when i run this on a sample of about 300 ways it looks like an infinite loop.

Does anybody have suggestions on where the actual problem is or have a better approach on finding subgraphs like the ones i am looking for?

For now i am working on the lines of Germany (a little bit filtered). But here i am looking for ways that match certain criteria (like "is wihtin a forest (polygon)") and now i want to build subgraphs out of them. Of course, i am well aware of this being ressource hungry. But in this scenario it would be ok to run for one or two hours.

Can @Mapper be a little bit more precise about that "smal set of subgraphs (Contraction hierarchies)"?

For each good_way that is not in another subgraph already i look recursively for intersecting private.good_ways and than update their subgraph_id in private.good_ways with a common id i retrieve from the sequence private.subgraph_seq. The table good_ways has the osm_id of the ways that intersect forrests and the like (that match criteria for good ways). The way of a good_way is retrieved via the join with planet_oms_line.
I hope you see what i try here. The problem seems to be, that the recursion is very very slow when the subgraph contains a lot of lines. But at least it should finish at some point, but it does not seem that way.


texts: ArcGIS Pro SDK set rotation / angle field in renderer

texts:
I'm in the process of converting an Addin written in ArcObjects to an ArcGIS Pro Addin. I need to rotation everything in a point layer using the values stored in the ANGlE field. In ArcObjects I was able to define a rotation renderer like this then apply to the entire point layer,

pRotationRenderer.RotationField = "ANGLE"    ' ANGLE is the name of a field in the feature class
pRotationRenderer.RotationType = esriSymbolRotationType.esriRotateSymbolGeographic

How do I achieve this using the ArcGIS Pro SDK?


texts: Generate grid for specific scale in QGIS

texts:
I have a rather large map that I want to break down into smaller parts, or sheets. So I'm going to use the 'atlas' feature of QGIS Print Layouts. My layers in QGIS have a CRS of GDA94 (degrees), where as my map layout will be projected to the specific UTM zone for that area (meters).

What I want to do is create a coverage layer containing a grid of polygons. This will essentially become the extents for each individual map sheet (as well as have other attributes such as map name, region etc). The issue I'm having is trying trying to determine the size of my polygons so that I can get it to be a specific size to match the size of the 'Map' object on my print layout at a specific scale, if that makes sense?

For example, my map sheets are standard A2 page size, which is 594mm x 420mm (portrait). Accounting for margins, let's say that the 'map' object on my print layout is 580mm x 360mm. What I need to do is calculate a grid of polygons where each polygon will be 580x360mm at scale exactly 1:50000 for that particular UTM zone. The scales are important so that I can utilise a grid where 2cm = 1km.

So far, when I try to generate the polygon grid, I only have the option of setting the size in degrees (as this is the unit of the layer CRS). I can get it close to 1.50000, but I dont want close enough, as this will not be precise enough.

Is there a way to calculate this or is my understanding of what I want completely off the path and I should be looking at a different solution?


texts: I am using Grass 7.2.1 as a package in QGIS(2.18.13). I want to use the r.skyview addons. I installed the library and tried to run it from Python in GRASS. When I run the code in the python shell, the dialog box (where I have to identify the input, output files and the parameters)  is not showing. as a result it gives error of missing input data. How to fix this.


texts:
I'm using the select by attributes to find unique values that I am copy pasting from a spreedsheet.

EX "jobnumber" = 1111116 OR "jobnumber" = 2222245 OR "jobnumber" = 5555655 .. etc

The problem is I have thousands of values and hit the character limit. I can do select by attribute again and just add to the selection, but this would require doing many, many times and would be pretty tedious.

I want to know if there is a better way of getting large numbers of unique values selected all at the same time, maybe using ArcMap's Python console, or making the spreedsheet into a table in ArcMap and relating it to the target layer.


texts:
I am trying to connect to a WMS server from this location through QGIS. I get the message

  The server you are trying to connect to does not seem to be a WMS
  server. Please check the URL

I am trying to connect by adding a layer from a WMS server.


texts: Why doesn't dissolve finish? It just stops without error message

texts:
I have a shapefile of the whole UK that contained Postcode area, district, and full postcode.

I was going to try and dissolve the regions into separate business areas by joining a lookup table that identified, for each postcode which business area it's assigned to. But I couldn't get the business area to appear in the dissolve popup window.

So instead I modified the DBF file, adding the new business area field. Now when I go into the dissolve function, the dissolve field does appear there but I can't get it to successfully finish. It just stops. No error messages. Nothing.

I think what I'm trying to do is simple but QGIS just doesn't seem to be happy.

This is the very first thing I'm trying to do with QGIS so please be patient with me.


texts: Iterating over point layer and select another layer's points in radius using PyQGIS?

texts:
I found this code snippet in Selecting features within certain buffer distance of selected feature using PyQGIS?, which is almost what I'm looking for, except I don't have polylines but points and I want to find another layers points in radius.

I tried to edit the code above, but I'm confused when to switch layers and how to select the points in radius properly, so that I can call my "do something"-function.

Here my code so far:

def selectPointsInRadius(bufDist):
  inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('points')[0]
  print ('points-layer found')
  #stores features in a dictionary for fast access
  featDict = {feature.id(): feature for (feature) in inputlayer.getFeatures()}
  #spatial index for faster spatial queries
  featIdx = QgsSpatialIndex(inputlayer.getFeatures())
  try:
    itlayer = QgsMapLayerRegistry.instance().mapLayersByName('itpoint')[0]
    print ('itpoint-layer found')
    itlayer.selectAll()
    for point in itlayer.selectedFeatures():
        xyList = point.geometry().asPoint()
        print(xyList)
        minx = min([xy[0] for xy in xyList])   #minimum x-coordinate
        maxx = max([xy[0] for xy in xyList])   #maximum x-coordinate
        miny = min([xy[1] for xy in xyList])   #minimum y-coordinate
        maxy = max([xy[1] for xy in xyList])   #maximum y-coordinate
        ids = featIdx.intersects(QgsRectangle(minx - bufDist, miny - bufDist, maxx + bufDist, maxy + bufDist))
        #lyr.selectByIds(ids)  #use this if you don't need total accuracy, much faster!!
        tids = []  #true ids of nearby features that we want to select
        for id in ids:
            if featDict[id].geometry().distance(line.geometry()) < bufDist:
                print ("found features in radius")
                tids.append(id)
        print ("selected features in radius")
        lyr.selectByIds(tids, 1)   #add to selection,
        #call doSomething-function here?
  except IndexError:
    print ('itpoint does not exist')


texts: one mxd per region or many mxd's per map

texts:
I am working somewhere where they insist on using one mxd per region, in this mxd there is every layer required to make maps (about 10 per region). I find it tedious because you have to turn layers on and off and recreate the legend, and I find that since each region has the same layers if someone has another mxd open I can't edit layers.

What is the best practice? the idea was that there were too many mxd's being created, but I always thought it was better to have one mxd per map.

we can't do data driven pages because each map has about 3 layouts and requires some fiddling that can't be automated.


texts: Seeking data source for worldwide cross-sections of topographic and bathymetric for woodworking art project

texts:
I am a decent R programmer and statistician, I am able to work with data sets pretty well.  I am less familiar with GIS tools, but if given a shapefile I can certainly plot political borders and differently-colored geographies with the mapping packages in the R language.
I have a pretty awesome woodworking project in mind, but I am not sure of a suitable data source, and I seek guidance about how to get topographic maps that do not end at sea level.  I just bought some brand new chisels and I'm going to buy a big rectangular slab of mahogany to start carving a map into.  title of the project: the world without water
I want to print a cross-sectional / overlay that would have been used to make these maps, except a worldwide version.  I apologize if that's the wrong terminology.
My woodworking idea is to relief carve -- Google image search -- a map of the world from the tallest mountains down to the Mariana trench.  that's 9,000m above sea level to 10,000m below sea level.  If I go by 500 meter intervals, that is exactly 39 cross-sections of the topography of earth that I would have to carve.  cross-section is the important term here, because for what I am trying to do, I need the entire earth's landmass at or below a certain height.
So for example:  cross section #1 would be landmasses that round up to 9,000m instead of rounding down to 8,500m.  this is exclusively mount Everest, since K2 rounds down to the cross-section map at 8,500m.  that means I should start carving my block of wood everywhere except the peak of mount Everest.
according to Wikipedia, the second, third, fourth, and fifth highest mountains are all closer to 8,500m than they are to 8,000m -- so they would all be left alone for the second iteration of the carving.
This process continues - where I carve one 500m layer of earth at a time - until I hit the bottom of the Mariana trench, 37 cross-sections later.
to re-phrase: cross-section map #1 should show all points on earth lower than 8,750m,  cross-section map #2 should show all points on earth lower than 8,250m, cross-section map #3 should show all points on earth lower than 7,750m, and this should keep going down until the deepest point.
I know nothing about what USGS data are available, and there's too much of it (all with fancy names) for me to discern what type of map makes the most sense for my idea.
I just want to be able to print these cross-sections out, so I can trace them and carve down layer-by-layer.  As someone with minimal knowledge of GIS nomenclature and available data, it's difficult to search around and know what to look for.
Where can I find data to create 39 cross-sectional layers of the earth's surface?


texts:
I have a problem with reading laz files that are stored at IBM cloud object storage.

I have built pywren-ibm library with all PDAL requirements - one of them with docker. I then deploy it to IBM cloud function as an action, where the error that appear is "Unable to open stream for 'Colorea.laz" with error 'No such file or directory.' How can I read the files with PDAL in cloud function?

Here is some of the code:

import pywren_ibm_cloud as pywren
import pdal
import json

def manip_data(bucket, key, data_stream):

        data = data_stream.read()
        cr_json ={
             "pipeline": [
              {
                  "type": "readers.las",
                  "filename": f"{key}"
              },
              {
                  "type":"filters.range",
                  "limits":"Classification[9:9]"
              }

          ]
        }

        consulta = json.dumps(cr_json, indent=4)
        pipeline = pdal.Pipeline(json.dumps(cr_json, indent=4))
        pipeline.validate()
        pipeline.loglevel = 8
        n_points = pipeline.execute()

bucketname = 'The bucket name'
pw = pywren.ibm_cf_executor(runtime='ammarokran/pywren-pdal:1.0')
pw.map(manip_data, bucketname, chunk_size=None)
print(pw.get_result())

The code is running from local pc with jupyter notebook.


texts: Calculate NDVI using ArcPy Spatial Analyst

texts:
My goal is to calculate NDVI using map algebra syntax.

import arcpy
from arcpy.sa import *

red_band_raster = (r"C:\B3.tif")
nir_band_raster = (r"C:\B4.tif")

#Calculating NDVI float raster
ndvi_raster = Divide(
    Float(Raster(nir_band_raster) - Raster(red_band_raster)), Float(Raster(nir_band_raster) + Raster(red_band_raster)))

#Saving raster:
ndvi_raster.save(r"C:\NDVI.tif")

Unfortunately, instead of varying NDVI values I get constant raster filled with "1.000000" and NoData


texts: Overriding precision control within Reclassify Tool of ArcGIS Desktop?

texts:
I have a coding scheme in which I am trying to apply using the Reclassify tool prior to combining the different raster datasets using raster calculator.  I have used the Reclassify tool for other datasets that do not exceed the precision limit; however, I cannot seem to reclassify a new value over 10,000,000,000.  The maximum value I will need to reclassify to is 100,000,000,000.

I have altered the raster data into 64-bit double precision but cannot seem to override the precision control within the Reclassify Tool (will not exceed 9).

Is there anyway to override this control using ArcMap 10.3.1 to accommodate these larger numbers without using ArcPy?

If not, are there any suggestions on how I can make the Value equal those higher values?

I have been working within ModelBuilder.


texts: How to add Map Legend, Scale Bar and Select tool to an ArcGIS Server Service?

texts: I have published a service from ArcMap to my Server Account but I can't figure out how to add the legend, scale bar and select tool to the service. According to the published service it was built using ArcGIS JavaScript API. Will I have to write some code in order to display what I want?


texts: Mask raster function in R not working to extract values in one raster based on conditions in another

texts:
require: raster package

I am having trouble masking values of one raster using values in another.  This is similar to this question:

https://stackoverflow.com/questions/51716315/how-to-subset-classify-raster-based-on-another-raster-grid-cells-values

but it doesn't have the moving window (focal function) as a component of the question.  I have two rasters with the same dimension, extent, projection.

Raster "MAT2resampled" is a raster with a wide range of temperature indices

Raster "wintPCP" is a second raster with pixel values of 1, O, and NA

The only pixels I want from MAT2resampled are those with the value 1 in wintPCP.

I used this code:
PthreshWint2 <- mask(MAT2resampled, wintPCP)

When I plot PthreshWint2, it doesn't appear as though the MAT2resampled raster has been masked at all.  I think I'm missing something simple.

I have posted the two rasters here:
https://drive.google.com/open?id=118kCu6QpoXSYmA17PpODRhHDJefHPkx1

https://drive.google.com/open?id=1N062TiHTHBqRYrclRW8O87WY-MX8taQg


texts: Inconsistent spatial comparison operations (ArcMap)

texts:
I'm trying to find geometries who's boundaries are coincident to a line. I've mostly been able to achieve this with IRelationalOperator.Touches in my code. Ideally my relation description of "FF*T*****" would work. My question comes from me playing around with the relation description.

ESRI defines Touches as:

  Two geometries are said to touch when the intersection of the
  geometries is non-empty, but the intersection of their interiors is
  empty.

Touch also has the spatial relation code "FF*FT****"

From that description I would expect the relation description of "F********" to return at minimum what Touches returns. However it does not.

In the following image I am comparing a green polygon and a blue polygon to a polyline. Touches will return both the green and blue polygons, but "F********" only returns the blue polygon. Why does the Touches function think there is no relationship between the interiors of the green polygon and the red line, but my custom relation description thinks there is?


texts: Can I prevent QGIS raster clipper from changing output cell/pixel size?

texts:
I have performed a DEM extraction using the clipper tool in QGIS 1.9.0 Master, using a mask layer, the original resolution of the DEM is 20x20 but in the resultant DEM the new resolution is 19.98x19.98. Both in meters.

Is it normal? How can I prevent QGIS clipper tool to change the pixel size?


texts:
I would like to use the freehand digitization tool for QField, but this does not work. I want to use it on a Huawei Y6 with an electric stylus pen. The edit mode is activated, a polygon layer is selected,
Consider mouse as a touchscreen device is issued.
When I touch the screen with the stylus, the icon for the freehand digitizing tool appears, but I cannot click it or anything like that. If I move the pen on the screen, the symbol disappears again.


texts: Exporting each layer in map to separate image using ArcPy?

texts:
I have a mxd file with multiple layers I have to export to PNG/JPG. I would like to create a script which turns on layers one at the time before they are exported. So far I have come up with the following:

import arcpy, string

mxd = arcpy.mapping.MapDocument("current")

#Read input parameters from script tool
LyrList = string.split.(arcpy.GetParameterAsText(0), ";")

PNGPath = arcpy.GetParameterAsText(1)

#Turn of all lyrs in list
for lyr in arcpy.mapping.ListLayers(mxd, LyrList):
    lyr.visible = False

#Turn on lyr and export as png one by one
for lyr in arcpy.mapping.ListLayers(mxd, LyrList):
    lyr.visible = True
    arcpy.mapping.ExportToPNG(mxd, PNGPath+"\\" + lyr.name + ".png")
    lyr.visible = False

del LyrList
del mxd

For the LyrList I use a dialoge box with the parameter data-type as Layer.
Does anyone have an idea how to get this script working? When I use the ListLayers wildcard the script does work, but since I have a lot of layerfiles which all hace a kind of unique name that does not do the job for me.


texts: Dissolve polygons from layer1 witch overlaps with polygons of layer2 using GRASS

texts:
I want to dissolve polygons from layer1 witch overlap with polygons of layer2 keeping attributes table of layer1 and using command line.

What I've done :

v.in.ogr -o input=original_plygons output=layer1
v.in.ogr -o input=plygons_to_add output=layer2
v.patch -a input=layer2 output=layer1
(v.dissolve input=layer1 output=dissolved column=??? )

the problem is to dissolve you should chose a column ! and after the patch (merge) the overlapping polygons of layer2 doesn't have a common same value to dissolve them correctly...

Can you please help me to solve this problem ?


texts:
The documentation for ESRI's File Geodatabase API (downloadable from http://www.esri.com/apps/products/download/) says regarding its MultiPartShapeBuffer Class:

  These functions provide access to the shape buffer. Consult the
  extended shapefile format document for the buffer layout.

I found the "ESRI Shapefile Technical Description" here: https://www.esri.com/library/whitepapers/pdfs/shapefile.pdf, but it doesn't describe curves and some of the other geometry that a file geodatabase can contain.

Where can I can I find the "extended shapefile format document" to which they are referring, or another source that contains byte by byte buffer layouts for these additional geometry types?


texts: Using ArcObjects Sample for Custom vertex editing commands?

texts:
I downloaded this sample :
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#/Sample_Custom_vertex_editing_commands/0001000002pz000000/

it's a Dynamic Link Library.
using Visual studio i generated the dll file from this project :

the dll file is VertexCommands_VB.dll.

i imported the DLL to my application and it's works fine BUT When i want to execute the application in other computer* it didn't work

How do I generate the Dll on this computer* to work with it?


texts: Pyqgis vector layer rollback and commitChanges

texts:
I have problems with the vector layer rollback and commitChanges function, it doesn't matter wich one i call the vector layer is is Changed.Here is my code i hope someone can help me. This is my firt post and im new with qgis.

from qgis.core import QgsApplication,QgsVectorLayer,QgsVectorDataProvider
from PyQt4.QtCore import QVariant

if name == "main":

QgsApplication.setPrefixPath('/usr', True)
QgsApplication.initQgis()

vLayer = QgsVectorLayer(
            "/home/diademuertos/scripts/qgis/data/point.shp",
            "vLayerTest",
            "ogr")

vLayer.startEditing()

caps = vLayer.dataProvider().capabilities()

if caps & QgsVectorDataProvider.ChangeAttributeValues:
    attrs = { 0: QVariant(2),1: QVariant("B"),2:QVariant("Dont!!!")}
    vLayer.dataProvider().changeAttributeValues({ 1 : attrs })
else:
    print "OOPS"

r = False
if r:
    print "commitChanges"
    vLayer.commitChanges()
else:
    print "rollback"
    vLayer.rollBack()

QgsApplication.exitQgis()


texts: How to get a Layer as ILayer (ArcObjects)

texts:
I had asked this question previously which I think will work.  But as I am new to this, I can't figure out how to get a Layer as an ILayer object to pass to the method that's in my previous questions answer.

I have been playing with the Mark Cederholm's 'Snippets.py' that uses comtypes to drive ArcObjects.  That works well, but I don't know enough to venture far from the examples given.  It first gets an IMxApplication object, and I can't figure out how to get a Layer file as an ILayer from this object so that I can use the script I was shown in my previous question.

In the end I want to be able to change the Display Expression of some Table Views, and be able to grab them by their name in the TOC.

Note: If you know how to do this with C# or VBA, I can figure it out in Python

Update: Here is the answer in Python which relies on Mark Cederholm's Snippets.py and comtypes:

from Snippets import *
GetDesktopModules()
import comtypes.gen.esriCarto as esriCarto
import comtypes.gen.esriArcMapUI as esriArcMapUI
InitStandalone()
pApp = GetApp()
pMxDoc = CType(pApp.Document, esriArcMapUI.IMxDocument)
pFeatureLayer = CType(pMxDoc.FocusMap.Layer(0), esriCarto.IFeatureLayer)


texts: Wrong SLD style generated by QGIS for Geoserver

texts:
can you tell me what is wrong with following SLD style? I've generated it by QGIS (2.18.3). In QGIS styles works fine, but when I upload it into Geoserver (2.10.1) and click valide I've had a lot of errors:

Legend of SLD does not works too:

Here is my generated SLD:
http://pastebin.com/GLp4xfx1


texts: Open more than one GeoTiff file with Rasterio and PyCharm at once

texts:
I have just started to deal with Rasterio and programming. And have the following problem that I have 20 GeoTiff files of which I need the extent.
But I am not able to open and read more than 1 GeoTiff at once.
Is it possible with rasterio to open multiple files at once, like the feature ListRasters from arcpy.
This is my code
import rasterio

extent = rasterio.open(r'C:\ms4w\Apache\cgi-bin\data\....tif')

bounds = extent.bounds
xmin = bounds.left
ymin = bounds.bottom
xmax = bounds.right
ymax = bounds.top

rectangle = "%s %s %s %s" % (xmin, ymin, xmax, ymax)

print(rectangle)

This is my output.
C:\Anaconda3\envs\BPLAN_851.py\python.exe C:/Users/.../Rasterio.py

460788.94715415215 5427005.507855503 461014.2159758464 5427270.910360379


texts: Create polygon of fixed area based on cost/friction surface

texts:
I am using ArcGIS Desktop 10.4 and have all extensions. I have a shapefile of 494 points. I also have a cost surface/friction surface raster. I want to create a polygon of 13.4 hectares around each point based on the friction surface. In other words I don't want a simple buffer around each point, The polygon should look irregular as it should follow the friction values, so for example if there is a steep hill to the east the friction value will be higher.

I also want to produce another polygon with the same points and friction surface that will outline one hours walk from the point.


texts:
I'm pretty new to Spatial ETL.

In accordance with this question I have experience the geokettle v3 but unfortunately i'm having a several problems working with it.

There is some issues that aren't being tracked since that geokettle development seems  knowing that:

  ..the future version of GeoKettle will be more pluggable with Kettle and
  will not be anymore a friendly and spatially enabled fork of Kettle.
  Hence, it will be possible to add spatial extensions provided by
  GeoKettle to any Kettle/PDI 5.x installation.

(GeoKettle_readme.txt)

Another issue that i've noticed is that Geokettle only leads with some geometry type like point, line or polygon and doesn't work with multipolygon or multiline.

So my issue is how can i make a Spatial ELT to migrate a PostGIS database with multilinestrings to another one PostGIS database?

I would like to know your experience about Spatial ETL? Anyone have already made a similar process (A spatial ETL with geometry type like MutliLineString) ?

Share your experience with Spatial extension for Talend or another Spatial ETL tool.


texts: Mapnik performance - large polyline dataset

texts:
I'm working with a very large set of polylines - contour data - to be rendered with Mapnik. (The contours are generated using gdal_contour from SRTM data.)

At present I'm storing them in PostGIS and Mapnik reads them from there. However, I'm not convinced this is the best use of space or CPU.

Would there be a performance hit if Mapnik read the shapefiles directly (with an .index file, of course)?

What would be the best way to minimise storage requirements for such a large dataset? Are there any other formats I should be considering?


texts:
I am new to QGIS. With the help of videos I am able to find distance between points of one layer to points of other layers. If I need to automate this process by selecting round objects as one layer and all rectangular as other layers. How would I do it?

Any guidelines how to achieve this would be highly useful.


texts: Remoteness Measurements or Accessibility Models in Rural Areas in Africa

texts:
I have a points shapefile (villages) that are located in remote areas in Africa. I need to give each village a measure of "remoteness." I dont just want to do a "as the crow flies" calculation from the village point to the nearest paved road. Can anyone give me specific instructions on how to do so? I'm using ArcGIS.

I know there are things like Grass R walk module and Naismith's Rule/Dijkstra's Shortest path Algorithm and using some combination of slope and topo rasters, but very specific instructions are much needed. I have a point shapefile with all the village points and a roads shapefile (but the villages are obviously very far away from the road network for the most part - at least the road network shapefile that is available to me)


texts: Trying to make a fish diversity map using ArcGIS 10.7.1

texts:
I'm relatively new using ArcGIS. I want to make a fish species richness map (or fish diversity map; essentially tells you how many fish species are in each grid cell). Things I did so far:

Obtained the fish data (they are hundreds of layers; each layer contains one fish species)
I merged all the fish layers (because I want to make a 10*10 km grid)
Made a 10*10 km grid (or fishnet)
Spatial joined the grid to the merged fish layer (problematic, I often end up getting null value in the attribute table and I'm not sure how do I count the number of species in each cell using field calculator or other functions; species description info in the attribute table are in text format).

Ideally, I'd like to make a map just like this (https://biodiversitymapping.org/index.php/usa-fish/).
How do I use the spatial join function in Arc to show the number of species per grid (or species richness per grid)? Also, is it possible to add in a new field in the attribute table which can count number of species (for instance, if there are 3 different species in one grid, the new field will show the value of 3, something like that)?
Again, my goal is to make a map that can show number of species per grid cell.


texts: Arcpy New Mosaic to New Raster script error

texts:
I wrote a script to Merge Raster files, but unfortunately scripts does not work:

import arcpy
import os

workspace = "D:\Kartes\Ortofoto_published\2400"
sr = "D:\Maris\Salni\CoordinateSystem\LKS 1992 Latvia TM.prj"

rasters = []
walk = arcpy.da.Walk(workspace, topdown=True, datatype="RasterDataset")
for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        rasters.append(os.path.join(dirpath, filename))

ras_list = ";".join(rasters)
arcpy.MosaicToNewRaster_management(ras_list,"D:\Kartes\Ortofoto_published", "newraster", sr, "8_BIT_UNSIGNED", "1000", "3", "LAST","FIRST")

Show error (ERROR 000582: Error occurred during execution)


texts: Do ArcGIS SpatialReference object factory codes correspond with EPSG numbers?

texts:
I've got a project where I'm reading in geometry objects from text files to put them in ArcSDE. These geometry objects have an associated EPSG code, and I'm looking for a way to ensure they're in the right projection when I've finished with them.

Having taken a look at the factory codes that I'm interested in, (GCS GDA94; EPSG:4283, and VicGrid 94; EPSG:3111) for a project I'm working on it appears that the EPSG codes correspond directly to the ESRI factory codes for the spatial reference (here and here). This is great, as it means I don't have to worry about navigating to projection files to create SpatialReference objects, for example in the code below:

import arcpy

# Create a spatial reference object using a factory code
#
sr = arcpy.SpatialReference()
sr.factoryCode = int(epsg_code.lstrip("EPSG:"))
sr.create()

Unfortunately I can't find any reference in the ESRI documentation to say whether this holds true all the time or not - I'm guess it does, but I'd rather know for sure. Does anyone know where I'd be able to find out?


texts:
I have Yakkety installed on my PC.

I have installed Quantum GIS 2.18 using the installation process described in the link https://www.qgis.org/en/site/forusers/alldownloads.html#debian-ubuntu. I used the debian method.

I am trying to use the OpenLayers Plugin, which is trusted by the QGIS S/W.

While enabling the plugin, it screams that it requires QtWebKit, which although installed on the PC, is not found by QGIS. Apparently QtWebKit (as required by QGIS) is for QT4, but Yakkety is on Qt5.

Although this issue is reported as a bug here (https://hub.qgis.org/issues/15128), what is the solution?


texts: Extract values of pixels from each band of a multi-band raster with a shapefile in Python

texts:
I have downloaded multiple products from EO-1 satellite and I have created digital rasters out of them. Each of them contains 10 or 242 bands (those are hyperspectral). I have also created a shapefile for each of them with multiple points. Now I would like to create a table/array with Python containing values of a pixel from each band with coordinates from a shapefile. I would like to have all bands in separate columns and all points in separate rows.
I have managed to create such data frame in R, and I tried to do so with GDAL and numpy but I am far less familiar with Python.


texts: QGIS: Mergin two layers always creates new attibute field

texts:
I want to merge several layers into one. Each layer has a specific Attribute field (type: GG, DO, eg.). In the merged layer QGIS produced a new Attribute field for each type.

I want QGIS to but all the values in one Attribute field called "TYPE". All the values should be listed in that field.


texts: How to write tif file in java from another tif file using GeoTIFF?

texts:
I am looking for reading tiff file then processing its pixels and write it back again as a tiff file. Is it possilbe, if yes please provide some 5-10 line sample.

ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY
        .createValue();
policy.setValue(OverviewPolicy.IGNORE);

// this will basically read 4 tiles worth of data at once from the disk...
ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();
//gridsize.setValue(512 * 4 + "," + 512);

// Setting read type: use JAI ImageRead (true) or ImageReaders read methods (false)
ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
useJaiRead.setValue(true);

//reader.read(new GeneralParameterValue[] { policy, gridsize, useJaiRead });
GridCoverage2D image
        = new GeoTiffReader(f).read(new GeneralParameterValue[]{policy, gridsize, useJaiRead});
Rectangle2D bounds2D = image.getEnvelope2D().getBounds2D();
bounds2D.getCenterX();
// calculate zoom level for the image
GridGeometry2D geometry = image.getGridGeometry();

BufferedImage img = ImageIO.read(f);

Here I got bufferedimage which i will process for some functions then

I would like to write this tiff again as a new file and I will process it further.


texts: How to create GeoTiff from image and inertial data with GDAL

texts:
I am trying to position photos that have seconds of week (utc) in their name with a given inertial data from a file, the file has the following columns:

UTCTime   Easting   Northing     H-MSL   Roll   Pitch   Heading

Matching the photo UTC time with the record of the inertial file I am trying to build the geotransform, the problem is that it looks like there is a close relation with the degree of the rotation and the pixel size, when I change one of the two, the other changes. For example, with pixel size as 1m photos are correctly rotated and positioned, but below 1m pixel size, the rotation is wrong.

There is barely any information about the world file:

The world file is an ASCII text file associated with an image and contains the following lines:

Line 1: x-dimension of a pixel in map units (A)
Line 2: rotation parameter (B)
Line 3: rotation parameter (C)
Line 4: NEGATIVE of y-dimension of a pixel in map units (D)
Line 5: x-coordinate of center of upper left pixel (E)
Line 6: y-coordinate of center of upper left pixel (F)

I am applying the following geotransform (Note, the coordinate is at the center of the photo):

(A) pixelSize * Cos(Heading)
(B) Sin(Heading)
(C) Sin(Heading)
(D) -pixelSize * Cos(Heading)
(E) Easting - ((photoWidth / 2 ) * pixelSize)
(F) Northing - ((photoHeight / 2 ) * pixelSize)

How do I compute the geotransform from Easting, Northing, and Heading?

Note: Heading is in degrees and I convert it to rad. PixelSize is user input.

When the pixel size is less than 1, the image is not properly rotated when I load it on QGIS.

Example:

PixelSize = 1

PixelSize = 0.1 (The case that mess the rotation pixSize<1)

Im just changing the pixel size in the geotranform and the rotation gets messed up.


texts: QGIS 2.14 GRASS geoalgorithms don't work?

texts:
I instaled qgis 2.14 and grass geoalgorithms dont work. Only 160 geoalgorithms and should be 168. MsBox: Missing dependensy, It seems that GRASS is not correctly installed and configured in your system. Please install it before running GRASS algorithms.
Also v.net geoalgorithms missing.

I reinstall QGIS twice.

I cannot find:
QGIS geoalgorithms:

Smooth: for smoothing line or polygon layers.
Reverse line direction.

there are not on the list...


texts: Can you use arcpy.Copy_management() with in_memory tables?

texts:
I have a table in an ArcSDE GDB I'd like to step through using Python. For performance reasons I wanted to load it into memory like so:

arcpy.Copy_management('C:/GDB.sde/MyTable', 'in_memory\\MyNewTable')

When I tried this in IDLE I got:

===================  RESTART ===================

Running it in the Python window in ArcMap brought the whole program down. So my question is, can Copy_management actually be used to copy a table into memory? The Esri documentation is pretty vague on what kinds of data you can pass.


texts: OpenLayers 3 view.setRotation shortest rotation change direction

texts:
I am experimenting with OpenLayers 3 for the first time and am looking to follow the user's position and heading.

I'm following their heading with

deviceOrientation.on('change:heading', function(event) {

    var heading = event.target.getHeading();
    view.setRotation(-heading);
});

This seems to work reasonably well but every few seconds the view will do almost 360 degree spin when only moving slightly. I believe this is because the setRotation is going "the long way round" to the new heading value. Is there a way to ensure it rotates the direction which involves the shortest movement to the new heading?