texts: Metadata property 'system:time_start' not found after scaling MODIS NPP data?

texts:
I'm trying to re-scale an image collection using the solution provided here and then to reduce it using the timestamp as per this example.
The script fails with the error:

  NPP trend: Layer error: ImageCollection.reduce: Error in map(ID=2004_01_01):
  Image.metadata: Metadata property 'system:time_start' not found.

Can anyone tell me where I'm going wrong?

My script:

  var addTime = function(image) {
  return image.addBands(image.metadata('system:time_start')
    .divide(1000 * 60 * 60 * 24 * 365));
};

var scaleNPP = function(image) {
  return image.expression('float(b("Npp")/10000)')
};

var collection = ee.ImageCollection("MODIS/006/MOD17A3H")
  .filterDate('2004-01-01', '2010-07-01')
  .filterBounds(NNam)
  .map(scaleNPP);

var trend = collection.select(['system:time_start', 'Npp'])
  .reduce(ee.Reducer.linearFit());


texts: Break LINESTRING values in postgis into mutiple LINESTRINGs

texts:
Figure A below is the visual representation of what i have in my road data set. Line 1, Line 2, Line 3 in illustration are each represented as 1 row postgis DB inside a geometry column, some of these geometry fields are LINESTRING some of them are MULTILINESTRING. So i totally have 3 rows in database for Figure A.

What i would like to do is to convert this data set to Figure B. In Figure B Line 1 and Line 2 are each broken in to 2 lines since they intersect with Line 3.

After this operation i would like have 5 rows represented as lines in DB instead of 3. I have no clue on how to achieve this.

Figure A

Line 1-----------------------------------
              |
              |
              | Line 2
Line 3-----------------------------------

Figure B

Line 1.1 -------- --------------------------- Line 1.2
                 |
                 |
                 | Line 3
Line 2.1 -------- --------------------------- Line 2.2

Of course this can be done easily after examining the data manually. Problem is i have 30K lines which need to be converted to more lines like this.


texts: Remote Sensing Indices for urban expansion

texts:
Apart from NDBI, could you name a few indices for urban expansion? I have done some search on google but haven't found anything that could help me. Perhaps,

What I'm interested is also: road/infrastructure, parks, parking lots, commercial areas  development in the urban environment.

Edit:

I need to analyze some orthophotos of a city and to do that I need to apply some indices to my analysis. As for software, I'm going to have at my disposal: ERDAS, ArcGIS, Open Source.


texts: In QGIS there is an option to create a custom form that should make it easier to create/update attribute entries. So my question here is how is the relation between the two? If I change something in the attribute table, does the value in the custom form also change and vice versa?


texts: In PostGIS, how to count nearby points efficiently?

texts:
I have one dataset A with around 6K point geometries and another dataset B with around 600K point geometries.

I am trying to come up with an efficient way to update dataset A with an integer representing the number of objects from dataset B that are within a certain distance (100 meters) from the point in A.

The strategy I came up with did not complete overnight, so I am thinking I could be doing this more efficiently. (Or maybe this is just that expensive an operation?).

I am running this on my X220 i5 laptop with SSD and 8GB RAM, Ubuntu 12.10, PostgreSQL 9.1, Postgis 2.0.

This is the function I use. nbi.testpoints is the small dataset, nbi.bridges is the bigger dataset. Both tables have spatial indexes on them.

CREATE OR REPLACE FUNCTION getnearbycounts() RETURNS VOID
AS
$$
DECLARE
  node nbi.testpoints%rowtype;
BEGIN
  FOR node IN SELECT * FROM nbi.testpoints
  LOOP
    UPDATE nbi.testpoints SET closenbicount = (SELECT count(1) FROM nbi.bridges WHERE
      ST_Transform(node.geom, 3785) <#> ST_Transform(nbi.bridges.geom, 3785) < 100)
      WHERE nbi.testpoints.id = node.id;
    RAISE NOTICE 'node % done', node.id;
  END LOOP;
END;
$$
LANGUAGE PLPGSQL;

I am doing the transform to a meter-based coordinate system to be able to give a distance threshold in meters, but replacing that with an approximation in degrees and doing away with the two ST_Transforms does not make a noticable difference. I see a performance of 1 point processed per 10-15 seconds.


texts: How do I edit info feature attributes in GeoServer with OpenLayers

texts:
Can you help me?  I have a very difficult problem when doing project...

I have used OpenLayers and GeoServer. When I click in the map to get feature info, I receive the info for the feature and region clicked.

Example: RegionName:Hanoi,Population:500...

I want to know how I can change the attribute value "RegionName" from Hanoi to HaThanh and save it in Geoserver. After I load the map I wish to receive the changed information.

Example: RegionName:HaThanh,Population:500...

Help me! Thanks.


texts: How to add a simple interface to mapserver CGI application?

texts:
I have a simple mapserver test page which shows information from shapefiles (cities, rivers, lakes, points of interest) etc. The result is a static PNG image. I want to add a simple interface so users would not have to modify the URL manually.

Options that I want to provide:

zoom in;
zoom out;
reset the map (back to default view with no zoom);
center map;
legend;
scalebar;
reference map;
show information about an object user clicked.

I know that map properties can be changed via URL but I am not sure about the way to do it. What I'm thinking is adding buttons in the template file and using javascript onClick events. When buttons are clicked, JavaScript redirects happen and additional parameters and parameter values are assigned to the URL.

I have read http://mapserver.org/cgi/ and this is what I have come up with.

In the example given below, button "Change map size" succesfully changes map dimensions to 800 by 400. I know that a full page reload happens, that is ok with me.

Is this the intended way of providing a simple GUI to a mapserver page? I want to stick to mapserver CGI, and not use PHP/Python mapserver API.

template.html:

<title>Mapserver Template HTML Document</title>

Template

<input type="button" onclick="zoomIn()" value="+" id="zoomIn" />
<input type="button" onclick="zoomOut()" value="-" id="zoomOut" />
<input type="button" onclick="reset()" value="Reset" id="reset" />
<input type="button" onclick="switchRegionVisibility()" value="Disable regions" id="switchRegionVisibility" />
<input type="button" onclick="changeMapSize()" value="Change map size" id="changeMapSize" />


texts: Display all mapfile layers with Openlayers

texts:
I'm using the following code to display my mapfile layers:

var wmsSource = new ol.source.ImageWMS({
            url: 'http://localhost/cgi-bin/mapserv.exe?MAP=/ms4w/apps/GEO/test.map',
            params: {'LAYERS': 'layer1,layer2...layerN'},
            serverType: 'mapserver',
            crossOrigin: 'anonymous'
});

Is there any simpler solution instead of typing all the layer names?

I mean something like this:
    params: {'LAYERS': all} ?


texts: ArcGIS sync replica not working with REST API

texts:
I have created a replica using CreateReplica REST API for offline use. Data format is SQLite. Now I am trying to sync replica.
Steps I have followed for sync:
 1. Uploaded my geodatabase file.
 2. Used upload id item in SyncReplica REST API.

Error:
Sync replica gives error as - can not import data.

I am not getting what is exact issue. How can I figure out the issue?


texts: Create vector layer from feature collection in OpenLayers

texts:
I have used a WPS process gs:UnionFeatureCollection to compute the geographical union of two collections of WFS vector features. The process returns a collection of features which I intend to display on the map, first, by constructing a layer and adding it to the layer switcher.

 wpsClient.execute({
    server:'local',
    process:'gs:UnionFeatureCollection',
    inputs:{
        first:features1,
        second:features2
    },
    success:function(output){
        var features = output;
        localStorage.setItem("Feature", JSON.stringify(JSON.decycle(features)));
        window.location.reload();
    }
})

This is the code to execute the process. I am saving the resulting feature collection to browser LocalStorage and refreshing the browser window. During the loading a new vector layer is being created and added to an array of layers that will be added to the map.

if(localStorage.getItem("Feature") != null){
    var features = JSON.retrocycle(JSON.parse(localStorage.getItem("Feature")));
    var vectorLayer = new OpenLayers.Layer.Vector("Union Layer");
    vectorLayer.addFeatures(features);
    layers.push(vectorLayer);
}

Here is my problem, the layer does appear on the layer switcher but not on the map even when enabled. Is defining features for the layer sufficient to display it or do I need to add it to the geoserver? (I'm pretty new to web mapping but not so much to javascript)


texts: Extent parameter PyQGIS gdal processing rasterize

texts:
I am having trouble with using gdal rasterize under windows in a PyQGIS script.
The command runs fine under Linux without specifying any RAST_EXT parameter. On a windows platform, however, the parameter apparently needs to be explicit (I get an error for the missing parameter value).
The problem is that I don't manage to pass this argument to the command. I get every time an error of wrong parameter value.
I have seen the posts on this topic, including:
How to pass an extent parameter to a Processing algorithm in PyQGIS?
How to make QGIS processing algorithms use default parameter?
Calling gdalogr:rasterize algorithm in python console in QGIS?
None of the answers are working for me.
My code looks like this.
#Import the grid
pointsLayer = QgsVectorLayer(nm, "region", "ogr")

#Define resolution in meters
resGrid = 2000

#Create the name of the output raster
outName = re.sub(".shp",".tif",nm)

extent = pointsLayer.extent()
extent = (extent.xMinimum(), extent.xMaximum(), extent.yMinimum(), extent.yMaximum())

#Save to a temporary file
tempout = str(os.path.join(tempfile.gettempdir(),"temp.shp"))
QgsVectorFileWriter.writeAsVectorFormat(pointsLayer, tempout, "utf-8", None, "ESRI Shapefile")

first = processing.runalg("gdalogr:rasterize", {"INPUT":tempout, "FIELD":"serp_id", "DIMENSIONS":1,
        "WIDTH": resGrid, "HEIGHT":resGrid, "TFW":False, "RTYPE":2,
        "NO_DATA":0, "RAST_EXT":'%f,%f,%f,%f' %extent, "BIGTIFF": 2, "OUTPUT":None})
first = first['OUTPUT']

I tried to provide the extent in several ways
"RAST_EXT": extent
"RAST_EXT": '%f,%f,%f,%f' %extent
"RAST_EXT": '%f, %f, %f, %f' %extent
"RAST_EXT": '%f %f %f %f' %extent

I also tried to switch the order of the extent values to match 'xMin, yMin, xMax, yMax', with again all of the combinations above.
To no avail.
As the name indicates pointsLayer is a points shapefile (grid of points), it is in UTM projection.
Again, this code with the same shapefile works great under linux without specifying a RAST_EXT parameter.


texts:
The way I see it, if two features to be labeled are too close to one another, then there are essentially three options. By default QGIS will omit one of the labels. If I enable the “show all labels” property, either for the layer in question or for the global labeling settings, then it will print overlapping labels.

Unless I have assigned data columns to store label positions, then I can resolve the conflict by moving one of the labels to a better position.

My problem is that I consider neither omitted labels nor overlapping labels an adequate solution. Is there a way to systematically iterate over all conflicting labels, giving me a chance to manually resolve those conflicts by moving one of the labels, or perhaps by hiding one of them?

As the placement seems to depend on scale and map extent, this might be specific to a given map in the print composer, I guess.


texts:
I need to be able to add bullet points to a text box so I can display a list on my map.  The work-around that I tried was to format the font differently and use a Wingdings character that looks like a bullet point.  This will work in the meantime, but instead of a top-to-bottom centered dot, the dot was situated more to the top.  It will work if need be, but I thought maybe someone else might have encountered the same problem.  Any suggestions would be helpful.

Using 10.0 desktop


texts: How to include external libraries with my QGIS plugin?

texts:
I'm developing a plugin which makes use of "networkx" library extensively. However, networkx is not part of Python libraries pre-installed in QGIS.

Is there a way to include this library with my plugin installation, so that my plugin can be portable?

Thanks!


texts: Comparing value with value from the next row

texts:
I have a table with two elements in this format:
a b  b a  b c  c b  c d ....f g g f

I have to analize if the second element is equal to the first element of the following row. If it's true, 1 is set in a new field. If not the marker is added by 1 +1. This is my code so far:
   with arcpy.da.SearchCursor(tab_selektiert, ["field1"],) as n1_cur:
    with arcpy.da.UpdateCursor(tab_selektiert, ["field2", "markerfield"]) as n_cur:
        filler_row = next(n1_cur)
        for n_row, n1_row in izip_longest(n_cur, n1_cur, fillvalue=filler_row):
            if n1_row[0] == n_row[0]:
                n_row1 = counter
                n_cur.updateRow(n_row)
            else:
                numcounter = + 1
                n_row1 = counter
                n_cur.updateRow(n_row)

del n_cur
del n1_cur

My problem is that beside the search cursor the n_cursor also skip one line. As result I only get b=b, c=c and missing a=b [..]
EditMy results should be something like this:  the field_2[0] and field_1[0](next row value) or field_1[0] and field_2[0](next row value) have to be equal (blue arrows). If so  a marker (started with 1) has to set in the markerfield. In the picuture  from line 6 to 7, the continuity makes a step forward and the marker has to update by 1 in line 7. Line 6 has '1' as marker even while the continuity breaks, because 'c' is in linie before.

Edit 2
counter = 1
numberlist = ["1", "2", "3", "4", "5", "6", "7", "8", "9", "0"]

with arcpy.da.UpdateCursor(tab_selektiert, (ufield1, ufield2)) as u_cur:
    with arcpy.da.SearchCursor(tab_selektiert, sfield1) as s_cur:
        next(s_cur)
        for sfield1, in s_cur:
            ufield1, ufield2 = next(u_cur)

        print("field 1  " + sfield1 + "---" + "  field 2   " + ufield1)
        print("field 2 + 1   " + str(int(ufield1[-1:]) + int(1)))
        print("fiel 1   " + str(sfield1[-1:]))
        print("---")

        if sfield1 == ufield1:
            ufield2 = counter
        elif sfield1[-1] in numberlist and ufield1[-1] in numberlist:
            if int(ufield1[-1:]) + int(1) == int(sfield1[-1]):
                ufield2 = counter
            elif ufield1[-1] == "9" and sfield1[-1] == "0":
                ufield2 = counter
        else:
            counter+=1
            ufield2 = counter

        u_cur.updateRow((ufield1, ufield2))

This gets me to the Point that all fields with a continuity get the marker 1. But a break in the continuity like 'd' to 'f' don't activate the else clause. This fields are just empty, also the last row and the following rows still gets the 1 as marker.


texts: Keeping symbol in Map but not showing it in Legend using ArcGIS API for JavaScript

texts:
I would like to customize the list displayed on the legend of the map but not delete it fro the Map.
Was able to get how to do it for desktop - but struggling to know how to do it in JS.
https://support.esri.com/en/technical-article/000016109
Am using ArcGIS v 4.15
this.renderer = {
  type: 'unique-value',
  field,
  legendOptions: {
    title
  },
  uniqueValueInfos: Object.entries(urlDict).map(([value, url]) => ({
    label: value,
    value,
    symbol: {
      type: 'picture-marker',
      url
    }
  }))
};


texts:
I have a next problem, when I convert my .shp file to .kml the polylines are all over the globe in Google Earth.

.SHP file is made in Autocad MAP 3D - export to shp. Note that MAP 3D is not creating .prj file so it cannot be opened in Google Earth pro.

I have made settings according to that question in the link below:

Problem when converting Shapefile to KML files in QGIS

Here are photos of input:

And here are photos of Google Earth output:

Any idea how to get proper output?

...I have opened .shp in QGis and saved it as ESRI .shp but Google Earth says it that it cannot be projected because coordinates are "out of area". I have set a WGS84 coordinate system in QGis.


texts:
I am using the OSGeo4W Shell with the ETOPO1, 1 arc-minute global bathymetric data (etopo1_ice_g_f4.flt).

I am cutting 256x256 pixel slippy map tiles, and am having a problem with getting smooth-appearing PNG files out of the GDAL process. For example, the tile at zoom level 10, x = 306, y = 389 represents a 24 arc-minute span, from 39.24 N to 39.64 N, and 72.42 W to 72.02 W.

So, to produce the GeoTIFF I apply:

gdal_translate -of GTiff -projwin -72.42 39.64 -72.02 39.24
-outsize 256 256 -r bilinear
F:\test\etopo1_ice_g_f4.flt F:\test\z10x306y389.tif

Then, to get the color relief tile, I apply:

gdaldem color-relief
F:\test\z10x306y389.tif F:\test\colors.txt F:\test\z10x306y389.png
-of PNG

The problem is that the resulting PNG file does not appear to have been properly interpolated. Instead, it is pixelated:

I have asked a similar question, thinking that the problem was because I was using a .NET wrapper for the GDAL libraries, but here I am using OSGeo4W Shell, and the results are the same: a pixelated image. So, this question is more about upsampling ETOPO1 data.

The pixelation remains the same regardless of the type of resampling method I choose (e.g., bilinear, cubic, lanczos, etc.), so I wonder if this is an issue with the ETOPO1 data source?

The intermediary TIFF file (z10x306y389.tif in my example) opens as "all black" in various image editors (e.g, Paint).

And the gdaldem hillshade option results in a simply terrible graphic:

Has anyone else had issues when upsampling ETOPO1 elevation data? Is there a way to check the "correctness" of the steps I am using, or the integrity of the intermediary GeoTIFF file?


texts: Combining attributes of point data within buffer using QGIS?

texts:
I have city data broken down into intersections and roads.  Around each intersection, I have created a buffer.  I also have parcel data, including the attribute "gross tax 2015" for each parcel.  I would like to assign each parcel intersecting a buffer to that intersection- and then analyze the total tax information for each intersection (there is more than one parcel in each buffer).  I have been able to create the buffer, and can visually see the parcels in each buffer, but can't seem to get an attribute table that shows to which intersection each parcel is assigned.  Furthermore, a parcel may reside in more than one buffer.  Is there any way to get the percent area of each parcel in each buffer?

I am using QGIS 2.8.6


texts: Proper projected coordinate system for Minnesota for measuring areas

texts:
I have seen in ArcGIS many different State Plane systems for Minnesota (location coordinates in WGS84: 44,07 lat/-94,61 long) within the Projected Coordinate system group:

NAD 1983 (CORS 96) (Meters)
NAD 1983 (CORS 96) (US Feet)
NAD 1983 (CORS 96) (Intl Feet)
NAD 1983 (Meters)
NAD 1983 (US Feet)
NAD 1983 (Intl Feet)
NAD 1983 HARN (Meters)
NAD 1983 HARN (US Feet)
NAD 1983 NSRS2007 (Meters)
NAD 1983 NSRS2007 (US Feet)

Which is the best choice for low area distortion?

Is that choice better than the UTM 15 N projection, for that purpose?


texts: python conditional statement to allow either two inputs or only one input going into a merge

texts:
I have a ModelBuilder process that I am converting to python for additional editing.  In the process I have both a point and line features(input variables) being buffered and then the resulting polygons merged into a single feature for additional processes.  But I want to make the model more versatile so it will run if only one of the inputs (point OR lines) is supplied.  Currently it crashes if only one input is supplied.  Is there a way to direct the model so that if only one input is provided, it skips the merge step (or successfully merges with null) and then goes to the subsequent steps.

Here's my code:

# Process: Buffer (6)
arcpy.Buffer_analysis(in_points_shp, point_buffer_shp, "buffer", "FULL", "ROUND", "NONE", "", "PLANAR")

# Process: Buffer (11)
arcpy.Buffer_analysis(in_LINE_shp, line_buffer_shp, "buffer", "FULL", "ROUND", "NONE", "", "PLANAR")

# Process: Merge (2)
arcpy.Merge_management("Z:\\supporting_docs\\point_buffer.shp;Z:\\supporting_docs\\line_buffer.shp", merge_output_shp,

...merge output then goes on to other processes.  My python skills are currently zero.


texts: ARCGIS-JS-API : featurelayer.queryFeatures() not returning results

texts:
I have a feature layer that is hosted on ArcServer. It has a field called document_num that is a string.

I have a table of which one of the columns is a link that lets the user "view on map". When this link is clicked it launches a function that is supposed to query the specific document by passing in the appropriate document_num..

The question is how should the where clause be formatted. I am having trouble getting the query to return any results. It is running though.

For example:
    If document_num = 100

function zoomToDoc(num){
    var query = new Query();
    var qString = "document_num" + num;
    query.where = qString;

    console.log(qTtring)

    (...)

    featureLayer.queryFeatures(query,......)
}

The console.log() in the above script produces "document_num=100". This runs successfully but returns no feature.

Should the where clause be structured such that quotes are wrapping the parameters? i.e. "document_num = '100'"   or    "'document_num'='100'"? I have tried all of the above but still return zero feature. And yes the number being supplied to the function is a valid entry.

Any thoughts as to where I am going wrong here?


texts:
I need to create tiles directly from c# code. The input file is in *.tif format. Is there a way to do this simply? I tried GDAL, but I found it complicated, I don't understand how to cut tiff into tiles (I'm novice in GIS :-(). Anyone can help me with logic how to create tiles for particular zoom level through GDAL? Probably, another way exists to do this task?

P.S. Unfortunately, using of gdal2tiles isn't suitable for me.


texts: Which Canadian province is better for studying GIS and finding professional opportunities: British Columbia or Nova Scotia?

texts:
I have a Business Administration degree and a Logistics MBA. My final project was a case study of eficiency improvements through ArcGIS routes in a local road transportation company. Suddenly I fell in love with GIS.

I actually live in Brazil, but I want to study GIS in Canada and, who knows, find a job opportunity to take my career to a next level. I'd like to stay in British Columbia or Nova Scotia because of the weather, more similar to Brazil. Both offer top rated programs.

In BC: BCIT, Okonagan, Selkirk College, SFU and VIU
In NS: NSCC (COGS)

I understand that BC offers more programs and more job opportunities than NS. The only reason I ask this question is because a cousin of mine is coincidentally immigrating to Halifax, NS. She also told me that Nova Scotia province has encouraging policies for international students and workers, specially in these areas.

So should I count on BC's hotter job market or on NS's conveniences?


texts: Create a turf polygon from extents to check intersection with other polygons

texts:
I am storing the bbox of layers as strings,
as the following string :

var layerExtent='[768039.2602094524, 5772524.37609651, 2313901.7202488566, 7435794.111581946]'

and I'm using this method to create a polygon from this bbox(extent) :

var bboxPolygon = turf.bboxPolygon(layerExtent);

but the polygon created has a weird coordinates so I can't use it to check intersection with other polygon with turf using:

var intersection = turf.intersect(bboxPolygon , item);

The structure of this polygon is as following:

so no intersection is detected, can you please give me ideas how make this work ?
can I create a proper polygon from extents so I could use it to check for intersection with other polygons ?


texts: Customized Albers Equal Area Conic projections are cosmetic only or do they affect area calculations too?

texts:
I will be doing area calculations for 100+ countries using country polygons taken from GADM, a database of global administrative areas. I am using ArcGIS. The GADM data is in WGS84 and is unprojected. I am doing some tests to gauge how area results differ whether I use customized projections vs. continental projections vs. national grids.

There does not appear to be a difference in calculated area whether I use Customized Albers Equal Area Conic Projection (I create country specific parallels and meridian) or Continental Albers Equal Area Conic Projection (built-in projections in ArcGIS). I thought a customized projection would lead to more accurate measurements but apparently the only difference is cosmetic - how the country is oriented on layout. In my limited test sample, all customized vs. continental results are the same except for Ghana. What might explain that? (Column F, Ghana)

I then compared the Albers resuls to results obtained using a country's national grid. There are very very minor differences vs. National Grids, always less than 1% except Brazil, where the difference is 1.12%.

What explains these minor difference and the Brazil difference in particular? Should the areas be the same?

Given that the observed differences are so small, I am inclined to do all calculations in Continental Albers Equal Area Conic rather than create customized parallels and meridians, or seek out national grids for all countries.

My core question is the whether Continental Albers Equal Area Conic Projections will always produce the same approximate areas as national grids (error <1%, column E)? Is Brazil the exception? Or is this a large country problem more generally? I did a few tests and do not want to test each and every country. I thought I could seek out the wisdom of this StackExchange.

I have scoured the Internet and can't find exact answers to these questions.


texts: Making point selection tool using ArcObjects with VBA?

texts:
I want to make a selection tool that I can use to select a feature point in ArcGIS for Desktop so I can combine it with a form that I just made to display the form after the feature selected.

I use VBA but, because I'm newbie, I don't really understand which interface that I should use.


texts:
I have heard that using ArcMap can replace some functions of ArcCatalog.

If it is true, how is the process done in ArcCatalog?

I don't have access to ArcCatalog all the time. I can't spend all my time in a computer lab just for ArcCatalog. I have access to ArcMap at all times.


texts: How to select a subset of a DEM in ArcScene?

texts:
I'm using ArcGIS 10.0 and trying to open a digital elevation model in ArcScene. I'm working on a project and got the DEM data from official sources which is quite huge. The resolution of the DEM is 1 meter. The total size is unknown but I think it's about 10km in diameter. The total file size is 1.5 GB.
Trying to add the DEM to ArcScene, using rightclick on Scene layers > Add Data and selecting the DEM, gives me the following message:

Could not add the specified data object to the scene. Try using a smaller subset.
Failed to open raster dataset.

Is the DEM too big? Why is that message? Could it be my total memory is not enough?
I'm running ArcScene 10.0 in a VirtualBox machine with Windows 7x86, with 2 processors 3.2 GHz and 4 GB memory allocated.
How to select a smaller subset? I neither see any option to do so nor find any subsets.


texts: Selected features disappear from map from ArcGIS API for JavaScript?

texts:
I am working with the Javascript API version 3.19. I have a feature layer used to display aerial images footprints. This feature layer is configured to fetch all the features to display at once based on an area of interest and a definition query using the mode FeatureLayer.MODE_SELECTION.

My issue is the following: the data is displayed properly at first, but if the user pans outside of the area of interest, all the selected features will disappear from the map, and will never show up again when I visit that part of the map again. From the debugger, I can see those features are still in the graphic collection associated with the feature layer, but they will not show up on the map. I have tried the following when the extent of the map changes:

Iterate through the features / graphics and call the method draw and / or show()

Call FeatureLayer.redraw()

Call FeatureLayer.refresh()

That still does not do the trick. Calling FeatureLayer.selectFeatures displays the data again but using that would make the code more complex than I suspect it needs to be. Does anyone have a suggestion?


texts: Organizing point labels using ArcGIS Desktop?

texts:
I am creating a map with about 50 points and want to label the points numerically (i.e. 1, 2, 3, 4...) and provide a table with the corresponding point numbers in order with additional columns with info related to each point.

My problem is I want the point numbering on the map to be grouped or organized so that the numbers are not randomly scattered all over the map (i.e. point 1 is near to point 2 etc) to assist in visual interpretation. My concern is if I just add a attribute field and asign numbers, on the map it will appear random.

Any thoughts on how to address this using ArcGIS Desktop 9.3.1?


texts: Is there a Java API for in-memory geo indexes?

texts:
My application is backed by PostGIS. There is a particular set of tables which have one of their columns being a geographic polygon. They are relatively small in size (few thousand rows) and can easily fit in memory. They do not change that frequently so I am thinking of caching them in memory and refreshing them every few hours in case there is a change.

The queries I will be getting is to find the polygon that contains a specific point. I might also need to get the closest polygon to the specific point (typically if it does not fall within any of the polygons in that specific table).

Is there a Java API that supports this functionality in memory? It will need to support the conventional geo indexes for fast lookup of these two queries I need to use. A lot of the Java GIS APIs I've seen seem to be backed by either Shape files or a database, so not really what I need, since I can already do that.

I was looking at this R-tree implementation: https://github.com/davidmoten/rtree
But for some reason it feels too low level.

On a separate question, I read this API is also possible to use, but it uses a quad-tree implementation: https://github.com/Esri/geometry-api-java

Apparently a quad-tree is more complicated (requires tile sizing configuration and similar fine-tuning parameters) and generally performs worse than an R-Tree. http://www.dpi.inpe.br/livros/bdados/artigos/oracle_r_tree.pdf

I was wondering if there is something which operates at a higher level, which supports the queries to get polygons that contain a point, and get closest polygon to a point, that is efficient.

UPDATE

Further search after a hint from one of the answers led me to GeoTools, which from the documentation seems quite comprehensive, which will probably be the first one to try out.


texts: Server (MapServer) does not support the gml coverage

texts:
I'm testing WCS (2.0.1) through OGC compliance test and some of the test cases are getting failed with an error:  The server does not support the gml coverage.

Test case id : cs2:get-kvp-core-req12

Assertion: req12: For all WCS request types, test the request version parameter has a fixed value of '2.0.1' by GetKVP.

Message d9183e2208_1:
    [MESSAGE] The server does not support the gml coverage.

Result: Failed

Any suggestions to solve this issue.


texts:
I am working with 2 adjacent 2014 LAS dataset provided in 2 projections:
However, when I type espg(las) my returned answer is  0
If I type summary(las) my coord.Ref.: NA

I wish to re-project both of these to NAD83/UTM18 and process them together through LiDR. I am unclear if this is possible to do in LiDR? How best to proceed? I have 1000 1km tiles to deal with.


texts: PyQgis convert screen to map coordinates

texts:
With PyQgis how can I convert absolute screen coordinates (say, 123,456) into map coordinates of the iface.mapCanvas?

It would be no problem if I had screen coordinates relative to the iface.mapCanvas.viewport(), since iface.mapCanvas.viewport().geometry() gives me sth like QRect(2,2,2198, 1090) and iface.mapCanvas().extent() gives me a QgsRectangle with xmin,xmax,ymin,ymax, so (leaving projection issues out of the scope) the conversion would be simple.

But since the absolute position of the mapCanvas may vary with toolbar positions or layer windows width, the question reduces to 'how to derive the absolute screen position of the map canvas?'.


texts:
I am attempting to select nulls from an in_memory feature layer using Arcpy's Select_analysis. Code example below, as well as, various iterations of attempts to alter the SQL query with no success.

sql_1 = '"{}_BBL" IS NOT NULL'.format(some_dataset_name)
print(sql_1)
arcpy.Select_analysis("Some_inmemory_layer", "Output_Name", sql_1)

The actual SQL query itself is as follows:

'"E_GIS_20190111_BBL" IS NOT NULL'

The empty output Feature Class yielded when encountering the invald expression error shows that the field name I am referencing is correct (field name matches above used in query in both attrib table and properties>fields tab).

I have tried the following fixes with no success:

'"E_GIS_20190111_BBL" <> NULL'

'"E_GIS_20190111_BBL" != NULL'

Error message received:

  arcgisscripting.ExecuteError: ERROR 000358: Invalid expression
  "E_GIS_20190111_BBL" IS NOT NULL
  Failed to execute (Select).


texts: Multiple attributes - is it possible to do variation by shape?

texts:
I am making a map of winter coat color for several species. I have two fields: species name and winter coat color. The latter can be either white or brown, and I would like to code, say, brown winter coat with circles and white winter coat with diamond, and then use colors to code species names. But Arcmap provides only two options: variation by color ramp and variation by size (see the picture).

Is there any way to do variation by shape in ArcMap?


texts: Swapping Lat and Lon coordinates of shapefile?

texts:
I have a shapefile that consists of one polygon whose boundary is defined by a large number of points.  Whoever originally made the file inverted the latitude and longitude coordinates, so the data displays at 38 degrees east of Greenwich and 122 degrees south of the equator, which is well south of the south pole, or off in space, or somewhere on the imaginary number line, or something.  How can I fix this?
I could make another shapefile and manually enter the correct coordinates, but this particular polygon is defined by several hundered vertices so I would rather not do it that way.  Is there a python script or GP tool that can help me?


texts:
I transform a csv file into a shapefile using the following script:

import shapefile as shp
import csv
import tkinter.filedialog
import urllib.request
import os
a = 'csvFiletoTransform.csv'

filename = 'C:/TEMP/ToscaAuxiliar.csv'
urllib.request.urlretrieve(a,filename)

out_file = tkinter.filedialog.asksaveasfilename()

nombre,codigo,y,x,time=[],[],[],[],[]

with open('C:/TEMP/xxxx.csv', 'r') as csvfile:
    r = csv.reader(csvfile, delimiter=',')
    for i,row in enumerate(r):
        if i > 0: #skip header
           nombre.append(row[0])
           codigo.append(row[1])
           y.append(float(row[2]))
           x.append(float(row[3]))
           time.append(str(row[4]))

w = shp.Writer(shp.POINT)
w.autoBalance = 1 #ensures gemoetry and attributes match
w.field('Longitud','F',10,8)
w.field('Latitud','F',10,8)
w.field('Date','D')
w.field('Date / Time','C',50)
w.field('ID','N')

for j,k in enumerate(x):
   w.point(k,y[j]) #write the geometry
   w.record(k,y[j],codigo[j], time[j], nombre[j])

os.remove('C:/TEMP/xxxx.csv')
 w.save(out_file)

it works, but the .prj file is lost.

It creates only the .dbf .shp & .shx files.

How can I do to create the prj file also ?


texts: Updating String title() for field in attribute table using ArcPy?

texts:
When using the title method, can it only be used on a string or can it run through an entire field and update Names to proper case in a table?

str.title()

I have tried the following and it did not work;

arcpy.title(table,'Field')


texts: Styling line layer based on two attributes in QGIS ? (1 discrete and 1 continuous)

texts:
I have found a really helpful tutorial to style points based on two attributes.

Is it possible to make style lines based on two attributes?

I have a shapefile with streets that have two attributes: traffic and lanes.

I can draw a thicker line based on the number of the lanes and color the lines from green to red with the traffic index, but I can't combine the two.


texts: Selecting parameters for visualization in Google Earth Engine

texts:
I don't know if I am the only one who struggles with this, but I seem to have a lot of problems selecting the best parameters to visualize images in Google Earth Engine. I am wondering if anyone has a methodology for doing this for multiband images (e.g. Landsat)?

I don't need to know what the best parameters are, but rather the best method for selecting parameters.

>     var ROI = ee.Geometry.Rectangle(-79.6179, 34.8945,-78.7994, 35.2923);
>
>     var landsat5 = "LANDSAT/LT05/C01/T1_SR";
>
>     var start = '1998-09-24';
>     var end = '1998-09-29';
>
>     var collection = ee.ImageCollection(landsat5)
>         .filter(ee.Filter.eq('WRS_PATH', 16))
>         .filter(ee.Filter.eq('WRS_ROW', 36))
>         .filterDate(start, end);
>
>     var image = ee.Image(collection.first()); // get first image
>
>     // Define the visualization parameters.
>     var visParams = {
>       bands: ['B3', 'B2', 'B1'],
>       min: 0,
>       max: 2500,
>       gamma: [0.95, 1.1, 1]
>     };
>
> Map.centerObject(ROI); Map.addLayer(image,visParams,
> 'image');

I know that I can change min, max, and gamma, but I always feel like I am just shooting in the dark.


texts: Importing from CSV and changing projection using ArcPy?

texts:
I'm attempting to convert CSV files to a group of shapefiles, then re-project them to a different coordinate system.

The first conversion works, but I'm having issues with the projection piece.

import arcpy
import os
arcpy.env.overwriteOutput = True
from os import listdir
from os.path import isfile, join, basename

csvpath = arcpy.GetParameterAsText(0)
outshapefolder = r'U:\skelley\Python\FWD_GPR_import\SampleCSVout2'

files = [join(csvpath,f) for f in listdir(csvpath) if isfile(join(csvpath, f))]

for f in files:
    arcpy.MakeTableView_management(in_table=f, out_view='tempview')
    arcpy.MakeXYEventLayer_management(table='tempview', in_x_field='Lat', in_y_field='Lon',
                                 out_layer='tempevent', spatial_reference=4326)
    arcpy.CopyFeatures_management(in_features='tempevent',
   out_feature_class=join(outshapefolder,basename(f).replace('.csv','.shp')))

#set new projection
outCS = arcpy.GetParameterAsText(2)
SR = arcpy.SpatialReference()  #an empty spatial reference object
SR.loadFromString(outCS)

inshape = r'U:\skelley\Python\FWD_GPR_import\SampleCSVout2'
outshape = arcpy.GetParameterAsText(1)

#loop through new shapefiles to project
for x in listdir(inshape):
    if x.endswith('.shp'):

        #Determine the new output feature class path and name
        infc = os.path.join(inshape, x)
        outfc = os.path.join(outshape, x)

        #Run project tool
        arcpy.Project_management(infc, outfc, outCS, 'NAD_1983_HARN_To_WGS_1984_2')


texts:
I have a large Raster with High resolution DEM of a large area, and I need to split it into smaller squares of 500mx500m, so that I can run other raster functions on it.

Basically I'm looking for the GDAL version of what I had asked here: Is there a Raster Equivalent of the Split tool?

How Do I Split a Raster into parts using GDAL?


texts:
I have a raster file (.jp2), I use pyQgis to set layer CRS (step1) and then save as a raster .tif with other CRS (step2). The result is not correct. It seems to be save as from original raster to CRS (step2) but I need the raster file save as from (step1) to CRS (step2).
Could you tell me how to resolve it?

Here is my code:

def getCustom(self, htd_output, f):
    htd_105 = "+proj=tmerc +lat_0=0 +lon_0=105 +k=0.9996 +x_0=500000 +y_0=0 +ellps=WGS84 +towgs84=191.90441429,39.30318279,111.45032835,0.00928836,-0.01975479,0.00427372,0.252906278 +units=m +no_defs"
    htd_106_nb = "+proj=tmerc +lat_0=0 +lon_0=106 +k=0.9999 +x_0=500000 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
    if htd_output == "UTM Zone 48":
        custom = htd_105
    elif htd_output == "VN-2000 Noi bo mui 3 KTT 106":
        custom = htd_106_nb
    res = QgsCoordinateReferenceSystem()
    res.createFromProj4(custom)
    return res
def save_file(self, type):
    layers = self.iface.legendInterface().layers()
    if type == 'tif':
        output_dir = self.dlg.lineEdit.text() + "/"

    htd = self.dlg.comboBox.currentText()
    htd_output = self.dlg.comboBox_output.currentText()
    # create directory if it doesn't exist
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
    for f in layers:
        if f.type() != 1:
            return False
        temp = self.getCustom(htd, f)
        if htd! = "System":
            f.setCrs(temp)
        if type == 'tif':
            extent = f.extent()
            width, height = f.width(), f.height()
            renderer = f.renderer()
            provider = f.dataProvider()

            pipe = QgsRasterPipe()
            pipe.set(provider.clone())
            pipe.set(renderer.clone())
            # pipe.set(f.crs())

            file_writer = QgsRasterFileWriter(output_dir + f.name() + ".tif")
            file_writer.writeRaster(pipe,
                                    width,
                                    height,
                                    extent, self.getCustom(htd_output, f))

And this's my UI in QGIS:


texts: Calculating incoming net radiation using Google Earth Engine?

texts:
I am a beginner with GEE and I am trying to calculate the incoming net radiation for that I need the square inverse of relative earth sun distance using the expression dr = 1 + 0.033*cos(J*2*pi/365) where J is the day of year.

Can someone help me with this code?

var sequence = ee.List.sequence(1, 365);
var dr2 =  sequence.divide(365).multiply(Math.PI).cos().add(1);
and alternatively I used
var addDate = function(image){
  var doy = image.date().getRelative('day', 'year');
  var doyBand = ee.Image.constant(doy).uint16().rename('doy');
  doyBand = doyBand.updateMask(image.select('B2').mask());

  return image.addBands(doyBand);
};
var image = ee.ImageCollection('LANDSAT/LT05/C01/T2')
  .filterDate('1987-01-01', '1990-05-01')
  .filterBounds(geometry);
var withDate = image.map(addDate);
print(withDate);
Map.addLayer(withdate);
var dr = withDate.divide(365).multiply(Math.PI).cos().add(1);


texts: Interpolation iteration using ArcGIS- arcpy

texts:
My current project requires interpolation of state level weather data. Approximately 5 years of daily data needs to be interpolated. I still need to figure out the best interpolation method. However, for this exercise I chose IDW. I found a script in ArcGIS forum and tried to use it. However, I am not sure why the script is not working. Any help would be greatly appreciated.

 #-------------------------------------------------------------------------------
    # Name:        Inverse Distance Method Interpolation
    # Purpose:
    #
    # Author:
    #
    # Created:     19/09/2013
    # Copyright:
    # Licence:
    #-------------------------------------------------------------------------------
    # ---------------------------------------------------------------------------
    # rainfall_idw2012.py
    # Created on: Tue Apr 08 2008 10:47:39 AM
    # ---------------------------------------------------------------------------

    # Import system modules
    import arcpy
    import sys
    import os
    import arcgisscripting
    import string

    arcpy.env.overwriteOutput =  True

    # Create the Geoprocessor object
    gp = arcgisscripting.create()

    # Check out any necessary licenses
    gp.checkoutextension("Spatial")

    # Load required toolboxes...
    gp.AddToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Toolboxes\Spatial Analyst Tools.tbx")

    #Output_variance_of_prediction_raster = ""
    modelPath = r"D:\Geodatabase\Rainfall\r2012"

    # Setting workspace
    gp.workspace = modelPath

    # Input Variables
    stations = r"D:\Geodatabase\MyGeodatabase.gdb\mr2012"
    studyArea = r"D:\Geodatabase\MyGeodatabase.gdb\boundary"

    try:
        # Get list of dau=ily events for looping
        daily_events = gp.ListFields(r"D:\Geodatabase\MyGeodatabase.gdb\mr2012", "r2012*")

        daily_events.reset()

        # Get the first reading and start loop

        dailyevent = daily_events.Next()

        while dailyevent:

            rain_amount = dailyevent.Name

            # Build output data layer name from the rainfall:

            r_level = modelPath + _rain_amount

            extract = modelPath + _rain_amount
            dailyevent = daily_events.Next()

            gp.AddMessage ("Current Operation: Preparing Interpolating using IDW")

            # Process: IDW interpolation

            gp.IDW_sa(stations, rain_amount, r_level, "1000", "2" "VARIABLE 12", "")

            # Process: Extract by Mask...
            gp.ExtractByMask_sa(r_level, studyArea, extract)

    except:
        gp.AddMessage("IDW failed " + gp.GetMessages())
        print "IDW failed ", gp.GetMessages()


texts: Adding processing buffer output to QgsMapLayerRegistry (argument 1 has unexpected type 'dict')

texts:
I'm writing a qgis plugin which needs a fixed-distance-buffer algorithm and must add the buffer result to the map layers of the software.

I have the following code:

# these are my import list
import os.path

import qgis.core
# Initialize Qt resources from file resources.py
import resources
import processing
from PyQt4 import QtCore, QtGui
from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication
from PyQt4.QtGui import QAction, QIcon

and the main part in run() method is:

    # region adding vector layers
    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        if layer.type() == qgis.core.QgsMapLayer.VectorLayer:
            layer_list.append(layer.name())
        else:
            break
    # region crime inputs and crime instance

    self.cInput = self.dlg.comboBox_input_Crime
    self.cInput.clear()
    self.cInput.addItems(layer_list)

after the if statement in run() method:

self.pic = self.cInput.currentIndex()
self.selectedCrime = layers[self.pic]
buff=  processing.runalg('qgis:fixeddistancebuffer', layers[self.pic], 100, 50, False, None)
qgis.core.QgsMapLayerRegistry.instance().addMapLayer(buff)

now when I run my plugin, after a progress bar appears, it shows an error of TypeError: QgsMapLayerRegistry.addMapLayer(QgsMapLayer, bool addToLegend=True): argument 1 has unexpected type 'dict'!

please explain what's the problem?


texts: Any example for a C#-WebGis development on the basis of DotSpatial?

texts:
I already used DotSpatial in a few projects. It is fine for reading/writing shapefiles, clipping etc..

According to the DotSpatial documentation, it is possible, to create WebGis projects, but I wonder why I cannot find any on the Web? Also, I wasn't able to find any tutorial for it. On GitHub there is a sample project called DemoWeb (https://github.com/DotSpatial/DotSpatial/tree/master/Source/Examples/DemoWEB) but I wasn't able to get it running.

Any hints?


texts: Label OSM polygons above a certain size?

texts:
I'm wanting to label polygons on a map that are above a certain area. These polygons are dervied from open streetmap and consist of either ways or multipart polygons ie If polygon is more than 160 sq m label it. I'm doing this using;

$area => 160

But nothing is being returned. The table has no specific area column but is comprised of polygons. On the fly transforms are enabled with all the source data in WGS 84 projected to a local national grid. The project CRS area is set to square metres as discussion in comments below.

What am I doing wrong?


texts: How to load a WMS layer in Google Maps API?

texts:
I wanted to expand my knowledge on the google maps API and get a hang of pulling in WMS layers. I have set up a simple map window and want to pull in a layer that contains electorate boundaries so they can be viewed.

I was looking through the documentation and tried my hand doing it as a GroundOverlay but cannot seem to get it working.

The WMS Url i'm trying to use is from:

http://data.gov.au/dataset/wa-state-electoral-boundaries-psma-administrative-boundaries/resource/b30c4ca8-f506-459c-b6d5-3df207e8666b

and my current code is:

function myMap() {
var mapCanvas = document.getElementById("map");
var mapOptions = {
        center: new google.maps.LatLng(-31.9505, 115.8605),
        zoom: 10,

    }
var map = new google.maps.Map(mapCanvas, mapOptions);

var overlay;
var imageBounds = {
    north: -13.67856638945,
    south: -35.1350604142437,
    east: 129.00192884357,
    west: 112.920934694525
    };

    overlay = new google.maps.GroundOverlay('http://data.gov.au/geoserver/wa-state-electoral-boundaries-psma-administrative-boundaries/wms?request=GetCapabilities',imageBounds);
    overlay.setMap(map);

    }

I guess I am just looking for a little prod in the right direction as to how to pull in a WMS with a URL like this. Am I looking at the wrong function entirely? I don't mind reading through the documentation and doing some self learning but just having trouble finding a starting point.


texts: Shifting ArcMap shapefile (lines) from point A to point B?

texts:
I would like to shift a shapefile (made of lines) from one point to another; not far away. This is really just for a cartographic adjustment, not a survey-level correction.

I was going to use Rubbersheeting but don't want to risk stretching the .shp; I only want to move it a bit.

I am just past a beginner level ArcGIS user.

I'm using ArcGIS 10.3.1.


texts: Dissolve Boundaries results in increased Shape_area than sum of parts using ArcGIS Pro

texts:
I have a series of overlapping/intersecting features from 3 years of data.
I buffered all feature types so they were polygons and then merged them. I made sure to map the fields. The buffers used were the actual values from the data and very nearly match those used in the excel calculations, though the stored value for the buffers was the max in ArcGIS Pro but we wanted to actually use the median in analysis. I am not the owner and cannot change the stored value.
Most features were very small and I don't think this explains the drastic increase.  It is usually a couple of meters max difference. When I dissolve the boundaries of this merged feature on species name, the shape area is far larger than the sum of the parts. Intuitively this number should be lower. For example: 3 polygons sit roughly on top of one another. The sum of the area for those 3 polygons is 100m.
So if I dissolve those 3 polygons into one polygon, shouldn't the area be about 1/3 of the previous sum?
In the original data (analyzed in excel and including sum of overlapping area) the total for this species is 10.4 acres. The sum for the buffered merged features is 42.6 acres, and the dissolved acres is 21.9 acres. See attached screenshots of 1) the merged data for a single species (yellow) and 2) the dissolved feature for that same species (red).


texts: Converting polygon feature class to WKT for ArcGIS

texts:
I am currently trying to export the geodata from a series of polygons in several .lyr files within my geodatabase in ArcGIS 10.2.
I want to export the series of polygons within my layer as WKT (well known text format). So that I may assign a geodata column in my CSV file for later importation and data management.

Does anyone know how to achieve this?


texts: How to extract cornerCoordinates with python and after convert to latlong?

texts: I need extract metadata from file .tiff but specifically I need get corner coordinates, and after convert this to lat-long.


texts: Disk access characteristics of accessing vs building an image pyramid.

texts: Hosting an extremely large (10-100x ram size) image pyramid via WMS has disk access characteristics closer to a random read, but generating the pyramid is a CPU bound sequential operation right?


texts: Measuring data load and redraw times in ArcMap?

texts:
Rephrasing of the question slightly to address some issues.

Is there are standard way of measuring load times for data, specifically sybolised data, within ArcMap?

In particular, we're comparing load times of different data formats (shapefile, SDE, file and personal geodatabases) across the network for the same (large) sets of data. We're also aiming to test speeds between ArcGIS versions.

At this point I've started hacking together a simple VB script to do some of the timing without having to sit there with a stopwatch as @matt-wilkie suggested. The script (thus far) can be found at http://textsnip.com/8912ac/vb. The script was written in ArcGIS 9.3 but works in 9.2 as well.

To use the script, copy the VB script to your mxd, and add two buttons, called "LoadDatasets" and "SymboliseDatasets". The LoadDatasets button allows the load of one or more feature classes or layers, and times the load. SymboliseDatasets checks the number of layers in the ActiveView, and if there are none, calls the load dialog (but doesn't time it). Once layers are added then the SymboliseDatasets button will symbolise all layers into 10 quantile groups based on their FIDs.

I have fixed the issue of timing the rendering by adding a DoEvents after the ActiveView.Refresh

In the meantime if anyone wants to pick this script up and modify it to make it more useful I'm happy to set this question to community wiki.


texts: Reduce overlapping / close polylines to one line

texts: I have a shapefile with multiple (often not perfectly overlapping lines). Attached is a picture of a typical situation found in the shapefile. I would like my output to have only one line (attributes are not important).


texts: WFS 1.1 - How to flip order of latitude and longitude in the response from Boundless Geoserver

texts: Please see the attached image. It shows a wfs 1.1 request and response. In the image I have indicated what I am trying to achieve in Red arrow. The coordinates array in the response has longitude first. But I need Latitude as the first item in the coordinates array in the response. Can anyone suggest me how to achieve that


texts: Creating a line along an existing polyline from a point

texts: I have a large data set of points, all of which lie on a river network and have x, y co-ordinates. I want to calculate the slope of the river channel at each point from the underlying DEM by taking the elevation at each point and the elevation of a point on the network 1km upstream. The best way I can think of doing this is by creating a line upstream from each point that stops when it reaches 1km away (obviously the line will not be straight in most instances), then create a point at the end of that line. I guess the flow direction or accumulation layer could be used to draw the line in the correct direction (upstream) from each point but I don't know where to start!


texts:
My goal is to create EPSG:4326 vector tiles (using Plate Caree) that can be served in Leaflet. I currently have a .shp file, which I converted to a geojson file using ogr2ogr. The main obstacle to begin with his seems to be how to create vector tiles from a shape file (or geojson).

Tippecanoe only produces Web Mercator tiles and the .mbtiles spec says reads "Only the Spherical Mercator projection is supported for presentation (tile display)". Therefore, I assume the preferred output is a MVT file.

The trouble is, I can't find any tool out there that can produce vector tiles in EPSG:4326. Does anybody know of a good way to produce a .mvt file (or a similar output format) in 4326 from a shape file?


texts: "The coordinates or measures are out of bounds" when converting labels to annotations

texts:
Am converting Labels to Annotations in ArcMap 10.1. However an error messages pops up with this message.

The coordinates or measures are out of bounds. Text = RUSHOZI

The coordinates or measures are out of bounds. Text = KIGAAGA

The coordinates or measures are out of bounds. Text = KAMUKUZI

Has anyone experienced this? How did you resolve this problem.

Please share the steps taken.


texts: Software that converts from legal description text?

texts:
Is there software that takes legal description text as input for conversion to some computer friendly GIS format? For example:

Those portions of the North Half of the Southeast Quarter and
the Northeast Quarter of the Southwest Quarter of Section 5,
T.13S., R.4E., S.B.B.&M., according to United States Government
Survey approved February 24, 1876 described as follows:

Beginning at a 1 1/2 inch iron pipe having coordinates
Y = 329,169.63, X = 1,895,486.55, set for the center of said
Section 5, thence along the East line of said Southwest
Quarter, S.4°11'18"W., 24.65 feet; thence S.56°02'45"E.,
...
All bearings and coordinates are referred to the California
Coordinate System, Zone 6.

More specifically, open source software.

Assuming that info that would be broad enough such as "S.B.B.&M." & "California
Coordinate System, Zone 6." can be addressed but realizing that a reference to something like a specific adjacent parcel recorded locally would require manual effort.


texts: Google Earth multiselect points to delete or drag & drop

texts: How can i select multiple items in Google's Earth sidebar "places" menu to delete or to move to a folder?


texts: SLD, combine variable and propertyIsEqualTo

texts:
I'm trying to make simple switch (if if else) statement in SLD with value coming from variable.

I'm calling Goeserver with parameter env=imageSize:A4.

Now based on that, I'd like to set cssProperty according to what value gets passed on. Pseudocode:

IF imageSize = A4
    cssProperty = 2
IF imageSize = A3
    cssProperty = 3
ELSE
    cssProperty = 1

I've come up with this:

<sld:CssParameter name="stroke-width">
  <ogc:Function name="if_then_else">
    <ogc:Function name="EqualTo">
      <ogc:PropertyName>pageSize</ogc:PropertyName>
      <ogc:Literal>A4</ogc:Literal>
    </ogc:Function>
    <ogc:Literal>2</ogc:Literal>
    <ogc:Literal>1</ogc:Literal>
    </ogc:Function>
</sld:CssParameter>

I can't get this to work, and after that I still need to implement second if statement.


texts:
Is it possible to cluster polygons in Leaflet?

In a similar way to the Marker Cluster plugin


texts:
I have a polygon layer and a point layer (with coordinates for each point).
The polygon layer got created using the MMQGIS plugin (Hexagonal Polygons).
I want to select all the points inside each polygon.

I worked with ArcGIS before but I want to change to open source.
Furthermore, I want to approach that through a Python script since I want to analyse the extracted data. I have QGIS 2.8 installed. I am a bit overwhelmed by all the GIS libraries there are for Python.

Can I approach that sort of selection with a GIS library that is already installed with QGIS? If so, could anyone pin point me into the right direction?


texts: Can MapboxGL use nearest neighbour interpolation for overlay image?

texts:
In HTML+CSS version of mapbox/leaflet I can set image-rendering option for .leaflet-image-layer with CSS to choose interpolation method. How can I make this in GL renderer?

Update: Solved

You can set nearest neighbour interpolation for all raster textures added to layers.

You need to change LINEAR parameter to NEAREST in 137 line of ImageSource

gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);


texts:
I'm looking for a mechanism/mechanisms where by I can pre-render my own map tiles and then server them up to the client through OGC WMTS service. There are numerous ways I have searched on the web which say how to get the pre-rendered map tiles and hand them over to a client, but none of them is talking about first how to pre-render the tiles themselves and also how to invalidate the old tiles and update the tile cache with new tiles.

And in my case WMTS service is a must and I don't want to take help of a WMS service to fetch a big map image and chop that up into small tiles and hand it over to the end client for users to see that.


texts: Is there a way to merge ImageCollections in Earth Engine without the "system:index" property being modified?

texts:
I was trying to merge two Sentinel-2 image collections (one based on date filtering and one based on image IDs), then pull in the matching Sentinel-2: Cloud Probability dataset to apply cloud masking to the merged collection.
However, it appears that performing a merge operation modifies the image "system:index" values (which are needed for the join with the Cloud probability dataset). The merge adds "1_" and "2_" prefixes to the "system:index" of the images in the merged collection. This stops any subsequent joining with another dataset from working.
The following is a simplified code of what I am trying to achieve. The problem is that the join with the Sentinel-2: Cloud Probability returns an empty set of images due to the prefixes added by the merge operation. If I simply join to one of the un-merged image collections, such as col1, the join works fine.
The "system:id" of the images are not modified by the merge operation and so I tried joining on the "system:id" instead of "system:index" but I received an error indicating:

FeatureCollection (Error)
Cannot filter on 'system:id'. Apply the filter to property 'system:index' instead.

Is there a way to merge image collections without the "system:index" properties being modified?
var testArea =  ee.Geometry.BBox(146.5, -18.3, 146.55, -18.2);

// Create two image collections. One based on image IDs, the
// second based on date filtering.
var IMAGE_IDS =  ["COPERNICUS/S2/20170104T002702_20170104T002921_T55KDV"];

var manualImages = IMAGE_IDS.map(function(id) {
    return ee.Image(id);
  });

var col1 = ee.ImageCollection(manualImages);

var col2 = ee.ImageCollection('COPERNICUS/S2')
    .filterBounds(testArea)
    .filterDate('2018-01-01', '2018-01-05');

// Verify the system:index prior to merge
var col1IDs = col1.aggregate_array('system:index');
print(col1IDs);
// ["20170104T002702_20170104T002921_T55KDV"]

var col2IDs = col2.aggregate_array('system:index');
print(col2IDs);
// ["20180104T002659_20180104T003132_T55KDV"]

// Merge the two collections into one collection
var mergedCollection = col1.merge(col2);

var mergedIDs = mergedCollection.aggregate_array('system:index');
print(mergedIDs);
// List (2 elements)
// 0: 1_20170104T002702_20170104T002921_T55KDV
// 1: 2_20180104T002659_20180104T003132_T55KDV
// The system:index has been modified to include 1_ and 2_ as prefixes

// Attempt to bind the S2_CLOUD_PROBABILITY dataset to these images
// Extract from https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_CLOUD_PROBABILITY
var s2Clouds = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
  .filterBounds(testArea);

// Join S2 with cloud probability dataset to add cloud mask.
var s2SrWithCloudMask = ee.Join.saveFirst('cloud_mask').apply({
  primary: mergedCollection,
  secondary: s2Clouds,
  condition:
      ee.Filter.equals({leftField: 'system:index', rightField: 'system:index'})
});

print(s2SrWithCloudMask)
// ImageCollection (0 elements)
// No match due to 1_ and 2_ prefixes from merge.```


texts: Making map in Robinson or Winkel Tripel projection in QGIS

texts:
I am trying to make a map like this:
https://docs.qgis.org/testing/en/docs/gentle_gis_introduction/coordinate_reference_systems.html#figure-robinson-projection

I have downloaded a shapefile with world borders from here:
http://thematicmapping.org/downloads/world_borders.php

This looks fine with the WGS 84 projection.

However, when I select the Winkel Tripel or the Robinson projection (both world and sphere), I get an artifact in Russia.

How can this be prevented?

I have located similar problems, but they are trying to project earth from space, which I am not, and therefore the solution does not seem to be applicable here.

Edit:
I have found one solution which is visually ok. I removed the part of the polygon that reached the edge of the map, on both sides. Mathematichally this is less than ideal of course.


texts: One element in attribute table not updating QGIS 2.18

texts:
When I use attribute table with option "Show Selected Features", and there is only one object selected, when I change values, they stay the same (table is still showing values before change) unless I unselect and select again the same object.
There is no problem when 2 features are selected (when I change value, the new value is showing immediately in table).

Is there any option to change that in QGIS 2.18?


texts: Getting rural/urban classification from Zip Codes

texts:
I have a dataset containing GPS coordinates and zip codes of locations from different regions in Spain, and I need to classify them as rural or urban based on either of these two information sets (and additionally having the city/town names based on these two).

I was wondering if anyone knows if there is any mapping (from zip codes or area names to rural/urban) I can find from a known, credible source, or if there is a way to classify these locations as rural or urban on my own if I have coordinates and/or zip codes?

Maybe if anyone can point me to a source that describes how this type of classification is typically done. For example, over here there is a mention of the classification that EU uses to classify as 'predominantly urban', 'predominantly rural' and 'intermediate', but is there a simpler method?


texts: How to Remove Decimal Points From Measured Grid Label C#

texts:
I'm creating a map layout with a measured grid using visual studio and ArcGIS Engine 10.4. I was able to change all the properties for the grid label (color, font and font size), except for the part after the decimal point. In other words, the changes I apply on the label using the code below are only applied on the part before the decimal point.

What I want to do is to hide the decimal part of the label. Anyone know how to access the decimal part and hide it or set the number of digits after the decimal to zero as in ArcMap for desktop?

 IGridLabel gridLabel = pMapGrid.LabelFormat;
 gridLabel.LabelOffset = 6;
 gridLabel.Color = color;

 gridLabel.LabelAlignment[esriGridAxisEnum.esriGridAxisRight] = false;
 gridLabel.LabelAlignment[esriGridAxisEnum.esriGridAxisLeft] = false;


texts: Applying the censusgeocode Package to an Entire Dataframe of Geocoded Data

texts:
I have a data frame of latitudes and longitudes. The goal is to apply the censusgeocode package to the entire data frame where our x is our latitude column and our y is our longitude column. After this, I am hoping to get a data frame with the geoid, state, county, tract, and block as columns.

Is it possible to do this so that I do not have to iterate over every single geolocation in a for loop? I could try the apply function, however, the apply function would not necessarily be able to be adaptable to this kind of problem. Below is my code for defining the dataframe and adding the solution to:

import pandas as pd
import censusgeocode as cg
data = pd.DataFrame({'latitude':[40.811765, 40.30019, 39.2464], 'longitude':[-81.528874, -81.83743, -83.60476]})

And what I have tried so far that takes too long to do this with thousands of locations:

geo_list = []
for f, g in zip(data['x'], data['y']):
    try:
        census = cg.coordinates(x=g, y=f)['2010 Census Blocks'][0]
        data = [census['GEOID'], census['STATE'], census['COUNTY'], census['TRACT'], census['BLOCK'], f, g]
        geo_list.append(data)
    except:
        pass


texts:
I am trying to speed up the process of loading labels in ArcMap. There are few layers containing labels that have a conditional form (not only in content but also text formatting). According to Esri guide:

  Label expressions slow drawing performance. In particular, avoid using
  label expression scripts to conditionally parse or format label
  strings or to concatenate the information in multiple fields into one
  label string when drawn. Instead, calculate a new field with exactly
  the information that is needed. The field can include formatting tags.

So instead of using label expressions I stored strings in Attribute Table.
I expected a lot faster loading, but labels are loading even a bit slower than before! What am I doing wrong?
Here is an example of label string:

<LIN leading = '-5'><FNT size='8' name='Arial narrow CE'><ITA>  So120-1.0</ITA></FNT> <FNT size='8' name='Arial narrow CE'><ITA>a </ITA></FNT><FNT size='8' name='Arial'>_____________</FNT></LIN> <LIN leading = '-5'><FNT size='8' name='Arial narrow CE'><ITA>1.63</ITA></FNT></LIN>


texts: Obtaining sink morphometry from local minima points

texts: I am working on LiDAR data from GIS and R and I want to extract landforms from the DEM. I have found a method to obtain local minima points of my LiDAR data. However, I would like to obtain the geometry of each sink. How can I do that? I have seen that Gradient Vector Flow algorithm (on matlab) and active contour can be used for it. Are there solutions to use it on GIS or R? Or another method?


texts: How to make transparent the background of topographic map in QGIS 1.8?

texts:
I have a raster topographic map, i.e. white background and black contour lines, houses, roads and so on. I want to make transparent the white background, so if I overlap the topographic map to an aerial photo, I see only contour lines, houses and other elements atop the photos.
In ArcGIS this operation is trivial. How can I do it on QGIS 1.8.0?
I explored many forums but I didn't find the answer. The solution is not the command "transparency" on the property window, because that command makes everything transparent.


texts: Rigorous geometrical models for InSAR imagery?

texts:
I'm looking for rigorous geometrical models for InSAR imagery. I'm familiar with the rigorous models for SAR images, but I don't know the model for InSAR. I searched through the internet, but I didn't find any clear-cut equations concerning this subject.

If anyone knows any books or lecture notes concerning rigorous geometrical models for InSAR imagery, please suggest me.


texts: Cannot edit QGIS project in QField (save to PostgreSQL db)

texts:
I struggle to edit my project in QField on a Samsung Galaxy Tab Active 2. The general idea would be to create and edit points in the field and save them in a PostgreSQL db (offline and online).

What I did so far is the following:

I created a portable project in QGIS (My project consists of a simple WMS base map and a "mask" that opens when creating a new point asking for the attributes of that point. This layer's source is the postgresql db. Editing and saving to the database works in the QGIS project.)
Exported the project via Qfield sync as a package for QField
Saved the package on the tablet under:
internaldata/Android/Data/ch.opemgis.qfield
Opened the project in QField on the tablet

-> The project will distplay perfectly but is not editable.

The Log-message says everywhere that the connection to database failed (PostGIS) and that editing cannot be started (QField).

To test this error, I set the layer in the QField sync project configuration to offline editing (so that it won't search for a db connection online). However, the log stays the same. This is why I believe that not being able to edit anything is an error that is based on some database connection mistake.


texts: Find maximum radius of circle that will fit within an irregular polygon?

texts:
I have a problem which I think could be handled by using the Zonal Geometry tool in the ArcGIS Spatial Analyst toolbox.  However I do not have a license for Spatial Analyst, so I am searching for an alternative; possibly using QGIS.

How do I find the maximum radius of a circle that will fit within an irregular polygon?

Note the polygon could be either a convex or concave hull (as shown below) and the solution must address both.

I tried Joseph's solution but unfortunately the result is not what I was looking for.

First, I do have very irregular polygons like this one:

If I follow Joseph's description the result looks like this:

This is for sure the result following that solution, but it is not the answer of my question.

Important for me is to answer the question how large the radius of a circle can be in maximum so that the circle is still completely inside the polygon, regardless of where the centre of the circle is.

For example, there is much more space in the north of the polygon, so that there can be placed a much larger circle than in the south of the polygon. But how large can this circle be?


texts: Exclude Values, including null, using Query Builder

texts:
I'm looking to exclude any streams that have values of "NULL" or "NP" using an SQL query in ArcGIS 10.2.

I've tried expressions like:

WTR_NM NOT (NULL, NP)

WTR_MN is the Field Name and I'm trying to remove the unnamed streams.


texts: Using SelectLayerByAttribute and Where Clause in ArcPy?

texts:
I'm trying to select the minimum value from a field to select and export that one feature to a temp shp so I can copy the field data from the point to a line. So far I have this:

import arcpy
from arcpy import env

# Set up the environment
env.workspace = r"C:\Users\BrettM\Desktop\OD_Python_Testing"
env.overwriteOutput = True

# Get Min
fcPoints = r"\OD_Shapefiles\One_Road_Points.shp"
temp = r"\OD_Shapefiles\intersectTempPoint.shp"

arcpy.Statistics_analysis(fcPoints, r"C:\Users\BrettM\Desktop\OD_Python_Testing\table",[["NEAR_DIST", "MIN"]])
nearIntersect = r"C:\Users\BrettM\Desktop\OD_Python_Testing\table"
minDistCursor = arcpy.da.SearchCursor("table", ["MIN_NEAR_DIST"])
for row in minDistCursor:
    intersectDist = row[0]

print intersectDist

arcpy.MakeFeatureLayer_management(fcPoints, "fcPoints_Lyr")
print "one"
arcpy.SelectLayerByAttribute_management("fcPoints_Lyr", "NEW_SELECTION", ' "NEAR_DIST" = intersectDist ')
print "two"
arcpy.CopyFeatures_management(fcPoints_Lyr, "\OD_Scratch\\temp")
print "three"

del row
del minDistCursor

All those print statements are just to see where I get in the script.  I get to print one.
I can't figure out how to get the where clause to work for the selectlayerbyattribute.  I keep getting an invalid expression error.  I may have the copyfeatures wrong as well but I haven't gotten to that part yet.

The overall goal is to cycle through a line file, select points that intersect that line and copy the data from the point that is nearest to a centerline to the line.  I'm using near to get the distance from the centerline for all the points and will work on adding the rest once I get this part to work.

There may be a much easier way to go about this but I'm very new to Python so I'm learning on the fly.


texts: Export from QGIS to Google Maps maintaining formatting

texts:
I'm totally new to QGIS so I'm not even sure this is possible.

I've created a bunch of polygons within QGIS and they look just like I want them (color, transparency, etc.).

What I'd like to do now is to export this somehow so that this same formatting is visible in a Google Map I'm creating as a layer. Is this possible? What format do I export?

I've tried exporting as geoJSON, but what appeared in Google Maps was not formatted. The boundaries of the polygons were there, but no colors, etc.


texts: Attributes missing from geodatabase in QGIS

texts:
I have been trying to load a gdb file into QGIS, but one of the layers do not show the attributes mentioned in the metadata. The data can be downloaded here. https://pubs.er.usgs.gov/publication/ofr20181088

The layer that doesn't work for me is called "Recent_2013_2017_Mining_Activity_Dataset". The metadata says there should be attributes such as the year of observation, ASM_val (which is the categorisation of mining activities) and others. However, all I am getting is a "count" attribute.

I really don't know what is causing this. Would you have any idea why these other attributes are not displaying in QGIS?


texts: Making correct shapefile for Mapbox Studio in QGIS?

texts:
I want to create an age building map of my home city. I made decision to follow this tutorial: https://www.wired.com/2013/10/how-to-map-building-ages-in-your-home-town/.

Building footprint shapefile for my city was not placed in open sources and I tried to use OpenStreetMap data. Via the OSM Downloader plugin I downloaded the multipolygons, filtered buildings, added required attributes and saved a shapefile.

I uploaded this file to this page to convert it to a Mapbox Studio tileset (a zip archive with files .shp, .shx, .dbf, .prj). I am seeing error Invalid zipfile (and no more information).

I reprojected layer to EPGS:3857, changed type of this layer to tile, converted geometry to the singlepart but this error still showing.

Trying to use other formats (i. e., kml) I get a tileset in Studio but my attributes were deleted.

How to fix this error? I use QGIS 3.8.0 (Zanzibar)


texts: Exporting Shapefile polygons from ArcGIS Pro to OpenStreetMap

texts: I recently worked on a project that digitized a lot of building footprints and swimming pools with high-resolution orthorectified imagery (~4 inches resolution). The polygons are all in shapefile format with a field indicating what the feature is. The polygons were created in ArcGIS Pro. I was wondering if there might be a way for this data to be contributed to OpenStreetMap since the buildings and pools digitized are not currently on OSM.


texts: How to put two points using query task please kindly solve the error in the following code

texts:
dojo.require("esri.map");
    dojo.require("esri.tasks.route");
    dojo.require("esri.layers.FeatureLayer");
    dojo.require("esri.geometry.Point");

    var map, routeTask, routeParams;
    var stopSymbol, routeSymbol, lastStop;
    function init() {
        map = new esri.Map("map", {
            basemap: "streets",
            center: [86.6611, 20.5848],
            zoom: 11
        });
        dojo.connect(map, "onClick", addStop);
        var featureLayer1 = new esri.layers.FeatureLayer("http://111.93.174.107:6080/arcgis/rest/services/network/Kendrapara/MapServer/1");
        map.addLayer(featureLayer1);
        var featureLayer = new esri.layers.FeatureLayer("http://111.93.174.107:6080/arcgis/rest/services/network/Kendrapara/MapServer/0");
        map.addLayer(featureLayer);
        //routeTask = new esri.tasks.RouteTask("http://192.168.0.5:6080/arcgis/rest/services/network/Network/NAServer/Route");
        routeTask = new esri.task.RouteTask("http://111.93.174.107:6080/arcgis/rest/services/network/KendraparaNetwork/NAServer/Route");
        //set up the route parameters
        //dojo.connect(routeTask, "onSolveComplete", showRoute);

        routeParams = new esri.tasks.RouteParameters();
        routeParams.stops = new esri.tasks.FeatureSet();
        routeParams.returnDirections = true;
        routeParams.outSpatialReference = { "wkid": 102100 };
        dojo.connect(routeTask, "onSolveComplete", showRoute);
        dojo.connect(routeTask, "onError", errorHandler);
        //define the symbology used to display the route
        stopSymbol = new esri.symbol.SimpleMarkerSymbol().setStyle(esri.symbol.SimpleMarkerSymbol.STYLE_CROSS).setSize(15);
        stopSymbol.outline.setWidth(4);
       routeSymbol = new esri.symbol.SimpleLineSymbol().setColor(new dojo.Color([0, 0, 255, 0.5])).setWidth(5);
    }

    //Adds a graphic when the user clicks the map. If 2 or more points exist, route is solved.

    function addStop(evt) {
        var stop;
        ////var x[]=[85.60, 20.65];
        //var myHome = { "lat": "85.60", "long": "20.65" };
        //var myHome1 = { "lat": "85.54", "long": "20.89" };
        ////            var x = "85.60, 20.65";
        //            //            var y = "85.54, 20.89"
        alert(xy);
        alert(x1y1);
        var start = new esri.geometry.Point(xy, new esri.SpatialReference({ wkid: 4326 }));
        var end = new esri.geometry.Point(x1y1, new esri.SpatialReference({ wkid: 4326 }));
        // alert(evt);
        for (var i = 0; i < 2; i++) {
            if (i == 0) {
                stop = map.graphics.add(new esri.Graphic(start, stopSymbol));
                alert("aaasaa");
                routeParams.stops.features.push(stop);
                alert("plz run yar");
            }
            else {
                stop = map.graphics.add(new esri.Graphic(end, stopSymbol));
                routeParams.stops.features.push(stop);
            }
        }
        if (routeParams.stops.features.length >= 2) {
            routeTask.solve(routeParams);
            lastStop = routeParams.stops.features.splice(0, 1)[0];
        }
    }

    //Adds the solved route to the map as a graphic
    function showRoute(solveResult) {
        map.graphics.add(solveResult.routeResults[0].route.setSymbol(routeSymbol));
    }

    //Displays any error returned by the Route Task
    function errorHandler(err) {
        alert("An error occurred\n" + err.message + "\n" + err.details.join("\n"));
        routeParams.stops.features.splice(0, 0, lastStop);
        map.graphics.remove(routeParams.stops.features.splice(1, 1)[0]);
    }


texts:
I want to find out the wkbType of a layer, this is my code:

print layer.wkbType()

Usually when it is a line geometry layer, it would give me 2. But instead it is giving me -2147483646.


texts:
I'm very new on GRASS and trying to extract the viewshed of a SRTM HGT file.

This is what I've done so far:
(why I'm using GRASS command line this way is not in question. Please answer this same way)

grass --exec r.in.srtm -1 --overwrite input=S23W044.hgt

grass --exec g.list type=raster

grass --exec g.region raster=S23W044 -p

grass --exec r.viewshed input=S23W044 output=S23W044-shed coordinates=-43.5,-22.5

But I can't find my output file and don't know it format.

GRASS GIS 7.5.svn


texts: Convert image from one EPSG defined CRS to another

texts:
I've got a lot of image files (.png or .jpg) that I want to convert from one EPSG defined CRS system to another, e.g. from EPSG:31255 to EPSG:25833, on a server within a Java app.

For each image I know the coordinates (x,y) of the lower left and upper right corner in that specific CRS system.

What I've found so far:

QGIS: Is able to do this but I'd have to do it by hand (answer)
GDAL: Is able to convert images but it uses C++ (which I can't use on the server)
Proj4js: Can only convert coordinates
Geotools: I've only found information about converting coordinates or full shapefiles so far.

Can GeoTools convert images too or what other Java/JavaScript library can?

Edit:

I'm currently using version 22.2 of GeoTools. The dependencies in the pom.xml file are:

junit
gt-shapefile
gt-swing
gt-geotiff
gt-image (not sure if I need that one)
gt-epsg-hsql

The imports in the java file are:

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.io.AbstractGridCoverage2DReader;
import org.geotools.coverage.grid.io.AbstractGridFormat;
import org.geotools.coverage.grid.io.GridFormatFinder;
import org.geotools.coverage.processing.Operations;
import org.geotools.gce.geotiff.GeoTiffFormat;
import org.geotools.gce.geotiff.GeoTiffWriter;
import org.geotools.referencing.CRS;
import org.geotools.util.factory.Hints;
import org.opengis.referencing.crs.CoordinateReferenceSystem;


texts: Multiple linear regression between grids: reviewing the assumptions

texts:
I am  performing a regression analysis with three raster stacks in R.
This worked so far.
Now I want to check the conditions of the model and I am not sure how exactly that should work in this case.
Since the model is calculated individually for each pixel, I cannot extract the residues for the entire investigation area and then check them.
Therefore I currently have to extract the data of the variables for a few points/coordinates on a trial basis and then create a linear model (lm) on which apply the tests (Shapiro-test, Breusch-Pagan-test, and VIF-test).

Does anyone know a more elegant solution?

NAO_stack <- stack("NAO_stack_fertig.tif")
AOD_stack <- stack("AOD_stack_fertig.tif")
CFC_stack <- stack("CFC_stack_fertig.tif")
NAO_AOD_stack <- stack(NAO_stack, AOD_stack)
NAO_CFC_stack <- stack(NAO_stack, CFC_stack)

NAO_AOD_CFC <- stack(NAO_stack, AOD_stack, CFC_stack)
NAO_CFC_AOD <- stack(NAO_stack, CFC_stack, AOD_stack)

fun1=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]);summary(m)$coefficients[,1] }}
Steigung_CFC <- calc(NAO_CFC_AOD, fun1)
names(Steigung_CFC) <- c("Intercept", "NAO_Steigung", "AOD_Steigung")
plot(Steigung_CFC[[3]])

fun2=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]);summary(m)$coefficients[,4] }}
p_value_CFC <- calc(NAO_CFC_AOD, fun2)
names(p_value_CFC) <- c("Intercept", "NAO_P_Value", "AOD_P_Value")
plot(p_value_CFC[[3]])

### pixel-by-pixel multiple linear regression with NAO and AOD / CFC as dependent variables

fun3=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]);summary(m)$adj.r.squared }}
r_squared_CFC <- calc(NAO_CFC_AOD, fun3)
plot(r_squared_CFC, main = "R²_CFC_multiple")

r_squared_AOD <- calc(NAO_AOD_CFC, fun3)
plot(r_squared_AOD, main = "R²_AOD_multiple")

### pixel-wise regression only with NAO as a dependent variable

fun4=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55]);summary(m)$adj.r.squared }}
r_squared_CFC_NAO <- calc(NAO_CFC_stack, fun4)
plot(r_squared_CFC_NAO, main = "R²_CFC")

r_squared_AOD_NAO <- calc(NAO_AOD_stack, fun4)
plot(r_squared_AOD_NAO, main = "R²_AOD")

r_squared_AOD_CFC <- calc(stack(CFC_stack, AOD_stack), fun4)
plot(r_squared_AOD_CFC, main = "R²_AOD_CFC")

### Extract data for a point to check the assumptions on the resulting LM model
XCoordinate = 8
YCoordinate = 61
points(XCoordinate, YCoordinate, pch = 15, cex = 0.5)
xy <- cbind(XCoordinate,YCoordinate)
xy
sp <- SpatialPoints(xy)
sp
data_CFC <- extract(CFC_stack, sp)
data_CFC
data_CFC <- as.vector(data_CFC)
data_AOD <- extract(AOD_stack, sp)
data_AOD
data_AOD <- as.vector(data_AOD)
data_NAO <- extract(NAO_stack, sp)
data_NAO
data_NAO <- as.vector(data_NAO)

CFC_lm <- lm(data_CFC ~ data_AOD + data_NAO)
summary(CFC_lm)
CFC_resid <- CFC_lm$residuals

AOD_lm <- lm(data_AOD ~ data_CFC + data_NAO)
summary(AOD_lm)
AOD_resid <- AOD_lm$residuals

shapiro.test(AOD_resid)
bptest(AOD_lm)
vif(AOD_lm)


texts: ArcGIS Raster To Multipoint Mac/OSX alternative

texts:
Is there an alternative to the ArcGIS Raster To Multipoint tool for Macs? The VIP thinning method is particularly of use. A comparative point cloud thinning/simplification tool to this would just as useful.

The objective is digital elevation model (DEM) (heightmap) to triangulated irregular network (TIN) (mesh) conversion on OSX.

Converting the DEM to a one-to-one (pixel to point) point cloud, to a TIN, and then using a mesh decimation algorithm is not a good option because of the sheer scale of the resulting point cloud and the time it would take to process it.


texts: Creating route between specific nodes based on sequence in QGIS

texts:
Right now I am working on a project related with the german railway network. I have a layer call "Stops_1_2" which contains the trip id, the stop sequence and the name of the station. I have another layer called "Railway_Germany" which I got from the Deutsche Bahn. I would like to create routes based on the common field trip_id and based on the network "Railway_Germany" so that I can get paths on the real network and not just straight lines.
I don't know how to adress this problem. I tried using shortes path ('Point to Point') but the problem is that I can only choose two points and I have a lot of data to do that manually.
This is the attribute table of my point layer....so what I need is for example create a for trip_id=10 following its stop_sequence, then create another route for trip_id=100 following the stop sequence and another route for trip_id=1001 following the stop sequence and so on that I cover all my trip_id.

This is my railway network and the stops

Is there a way to this on QGIS? Like I said I am new using QGIS and my programming skills are not good at all. (If there is maybe a solution with other software (like ArcMap) it would be helpful also)


texts: Not able to create a geoprocessing package in ArcMap using my own Network Dataset and associated data

texts:
I went through the following two examples to get an idea of how to create a custom, local Routing Task . . .

This is a walkthrough that takes you step by step through the creation of your first routing geoprocessing package:
http://resources.arcgis.com/en/help/...00000m5000000/
And here is the second part, that completes the picture and shows you how to run it from code:
http://resources.arcgis.com/en/help/...00000m6000000/

I am able to create the gpk using the above mentioned link instructions. However, when I attempted to create a gpk using my own data, I got an error message each time (General function failure ERROR 999999: Error executing function). When I researched this error, I didn't find anything helpful.
I am using a map document that contains a Network Dataset created from a Shapefile Feature class.
The settings that I use to create the package are . . .
Support ArcGIS Runtime enabled (other two options are disabled)
Summary and Tags detail filled in
Is there any way to find out more info on the error? I don't know why it's failing, but I keep changing little details and trying it again.

Thanks,
Renee


texts:
I have one gdb with lots of layer (~500) and to receive empty gdb I delete all records from layers. I Use deleterows / TruncateTable but for one layers take about 1-3 minutes.

Is there a faster way to get an empty gdb without records using python?


texts: I want to be able to run the NNJoin plugin via the QGIS python console to run a number of times for polygon-polygon and polygon-point (all shapefiles). I gather from reading a few other questions/answers around the subject that I first need to know the module? I would imagine there's a bunch of other code too? Can anyone help point me in the right direction?


texts: Exporting/Printing/Saving CMYK colored map to SVG, PDF, or PRN/EPS file

texts:
I am a beginner using QGIS 2.0.1 on Windows.

Following Data-defined Styles in QGIS
I fill CMYK values on my map layers, but I cannot Export/Print/Save CMYK to SVG, PDF, or PRN/EPS file. I see only RGB values on exported vectors.

How do I export CMYK colored map to prepress software like Illustrator?


texts: Converting WMS services into vector data using ArcGIS Desktop?

texts:
I have a WMS service which is distributed by the national organisation of my country. However, I want to process these data, but it's impossible due to WMS service is provided as a map.

Can I find any tool which allows me to convert the WMS service into vector data?

I have tried to export the WMS service into PNG image, specifying color mode as 1-bit Monochrome Threshold, in order to use it subsequently at arcScan tool.

However, the image which is extracted does not contain anything.

Is there any explanation?


texts: Inner shadow of a OpenLayers map container?

texts:
Not sure if this question should be here in the GIS section or if I should move over to a more common html/css place, but here we go:

I'd like to have an inner shadow of my OpenLayers map canvas so that it looks like that the map is underneath the rest of the web page. Example:

The square pattern is where the map will be later. Ideal would be to just use

box-shadow: inset 0 0 5px 5px black;

but the map tiles will be displayed on top of the frame which makes the shadow disappear.

Any ideas how to achieve this only with CSS? I really don't want to use lots of divs with PNG files as backgrounds, but I guess I have to?


texts: Organizing files for large Python script tool?

texts:
I am writing a lengthy ArcPy tool and am working with multiple python scripts, simply for organization's sake. Everything runs well when I run it from the command line, but now that I am starting to think about converting it to a script tool, I am getting a headache.

My previous structure:

parameters.py : lists all globals (parameters, functions, variables)
a_functions.py, b_functions.py, etc. : lists functions for use in main code, imports parameters with *
maincode.py : calls functions from *_functions.py, imports parameters with *

The structure is quite neat, and although it may not be ideal from a file structure standpoint, it works well enough (chime in if I am doing something horribly wrong).

Now that I need to load a file into the script tool that defines the parameters, I can't simply move all the parameters to the main code file, as I need to import them as globals to the other files, and importing with * will also import all the globals local to maincode.py.

I could just move everything into one file, but all the organization would be lost and it would be a big mess.

Does anyone have any suggestions who has built a lengthy arcpy tool with multiple files before?


texts:
I have Leaflet map with base layer (epsg4326) and wms layer of gps point (also epsg4326).
I've noticed when I zoom out , the gps point offset to the left.

Any idea why it happens ?

Here you can see that when I zoom out the symbol moves left each zoom.

This is the code:

function initLeaflet() {

    osm = L.tileLayer.wms(geoserver + "/OSM/wms", {
        layers: 'OSM',
        format: 'image/png',
        transparent: true
    });

    coords = L.tileLayer.wms(geoserver + "/coords/wms", {
        layers: 'Roham:DeviceInfo',
        format: 'image/png',
        transparent: true
    });

    var baseLayers = {
        "osm": osm
    };

    var overlays = {
        "coords": coords
    };

    map = new L.Map('map2d', {
        center: new L.LatLng(32.05522664582287, 35.05812072269409),
        zoom: 9,
        layers: [osm, coords],
        zoomControl: true
    });

    L.control.layers(baseLayers, overlays).addTo(map);
}


texts: Using DTM slope map to find shortest route between two locations using QGIS?

texts:
I've got two locations, A and B, and a DTM covering the region between A and B. Say I want to drive from A to B. I want the shortest route possible, but avoiding regions with slopes higher than, let's say 30 degrees.

Can QGIS calculate such a path?

If needed I can also script using python.


texts:
I would like make an animation in GeoServer using three shapefile.
I looked this page https://docs.geoserver.org/stable/en/user/services/wms/time.html#specifying-a-list-of-times but I didn't understand how I can create a "time specification".
How can I do this?


texts: Computing intersections for each geometry in a table

texts:
I have a table with many Linestrings and for each of these Linestrings I need to figure out which other Linestrings, from the same table, intersect with it.

If I recall correctly this gives me a Complexity of n^2. I tried doing this via a for loop in python and calling a SQL intersect() for each trajectory but it's taking forever: even after threading.

Is there a way of doing this computation in one SQL statement or perhaps using a different, less costly gis function call?

Here is the Python code (using psycopg2) that is called for each trajectory to compute other trajectories that intersect with it:

db_cursor.execute('''
                  SELECT
                      id
                  FROM trajectory
                  WHERE st_intersects((SELECT
                                            line_geometry
                                            FROM trajectory
                                            WHERE id = (%s)), trajectory.line_geometry);
                   ''', [trajectory[0]])

The database is PostgreSQL using PostGIS. Here are all the details for the trajectory table that I'm using: